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ABSTRACT 

Recovering the 3D structure of the scene from images yields useful information for tasks 
such as shape and scene recognition, object detection, or motion planning and object grasp- 
ing in robotics. 

In this thesis, we introduce a general machine learning approach called unsupervised 
CRF learning based on maximizing the conditional likelihood. We describe the application 
of our machine learning approach to computer vision systems that recover the 3-D scene 
geometry from images. We focus on recovering 3D geometry from single images, stereo 
pairs and video sequences. Building these systems requires algorithms for doing inference 
as well as learning the parameters of conditional Markov random fields (MRF). Unlike 
previous work, our system is trained unsupervisedly without using ground-truth labeled 
data. 

We employ a slanted-plane stereo vision model in which we use a fixed over-segmentation 
to segment the left image into coherent regions called superpixels. We then assign a dispar- 
ity plane for each superpixel. We formulate the problem of inferring plane parameters as 
an MRF labelling problem, which can be solved by an energy minimization method. The 
MRF is a graphical model in which superpixels define nodes and the adjacency between 
superpixels define edges. Our stereo energy function balances between a data matching 
term and a smoothness term. 

For systems with continuous valued variables, or discrete-valued variables with very 
large state space, it is impossible to directly use a standard discrete MRF inference tech- 
niques such as Loopy BR graph cuts or tree-reweighted message passing. For such systems, 
we propose to use a generic Particle-based Belief Propagation (PBP) algorithm closely re- 
lated to previous work, which we then formulate specifically for our MRF labeling prob- 
lems. Although we only describe a specific use of this generic PBP algorithm, we believe it 
can be used as an approximate inference scheme for a wide variety of problems that can be 
formulated by a probabilistic graphical model, especially those containing many random 
variables with very large or continuous domains. 

We demonstrate the use of our unsupervised CRF learning algorithm for a parameter- 
ized slanted-plane stereo vision model involving shape from texture cues. This unsuper- 
vised learning algorithm implicitly trains shape from texture monocular surface orientation 

iii 
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cues. We exhibit that training monocular cues from stereo pair data improves stereo depth 
estimation. Our stereo model with texture cues, only by unsupervised training, outperform 
the results in related work on the same stereo dataset. 

Our unsupervised learning method is also implemented for the monocular depth esti- 
mation (MDE) problem. The MDE model, learned using stereo pairs only, demonstrates 
a modest improvement after a few training steps, and achieve performance comparable to 
previous work on the same dataset. The use of MDE in combination with the dense stereo 
model also introduces a small boost in depth estimation over the initial stereo model. 

In this thesis, we also address the use of stereo video sequences. We formulate struc- 
ture and motion estimation as an energy minimization problem, in which the model is an 
extension of our slanted-plane stereo vision model that also handles surface velocity. Sur- 
face estimation is done using our own slanted-plane stereo algorithm. Velocity estimation 
is achieved by solving an MRF labeling problem using Loopy BR Performance analysis 
is done using our novel evaluation metrics based on the notion of view prediction error. 
Experiments on road-driving stereo sequences show encouraging results. 



ACKNOWLEDGMENTS 

It has been an honor for me to be advised by Prof. David McAllester, who has been giving 
me his outstanding mentorship, his patience and unconditional support from the first day to 
the last day of my graduate program at TTI-C. Through him, I also got to understand more 
about how to think and act as a true scientist. I would never have come this far without his 
guidance and encouragement. 

Many thanks go to my thesis committee members: Greg Shakhnarovich, Raquel Urta- 
sun, Ronen Basri for their valuable and constructive advices and feedbacks. 

During the last few years, I have received a lot of help and support from people in TTI- 
C, CS Department and elsewhere that I want to give thanks to: Pedro Felzenszwalb for his 
inspiration that fired up my passion in Computer Vision; Alex Ihler, Ronen Basri, Deva 
Ramanan, Xiaofeng Ren, Cristian Sminchisescu, Yali Amit, John Langford for helping 
me find answers to lots of my questions and puzzles; Wonseok Chae, Andy Cotter and 
Karthik Sridharan for being my friend and making me feel less alone as a TTI-C student; 
Gary Hamburg, Carole Flemming, Christina Novak, Adam Bohlander, Julia MacGlashan, 
Dawn Gardner and others, who have created a family-like environment and atmosphere 
at TTI-C, where I consider my second home away from home during the last few years; 
Nathan Srebro and Umut Acar who gave me advices as Academic advisors; Seiichi Mita 
and TTI-J Nagoya for welcoming me and giving an warm experience in Nagoya; Mr. Kien 
Pham and Adam Kalai, through whom I got to know about TTI-C and then became one of 
the first students here; Robinson, KC Lee, Katharine, Paul King, Dan Prochaska for their 
support and guidance during my summer internship. Special thanks go to my friends and 
my family. Last but not least, I want to thank my parents and my beloved Nhung. 

This thesis was partly funded by a grant from the Vietnam Education Foundation (VEF). 



TABLE OF CONTENTS 

ABSTRACT ED 

ACKNOWLEDGMENTS m 

LIST OF FIGURES EH 

LIST OF TABLES ED 

Chapter 

1 INTRODUCTION [Tj 

1 . 1 Markov Random Fields and Related Formalisms |2] 

1.1.1 MRF |2] 

1.1.2 CRF [3] 

1.1.3 Hidden CRF H 

1.1.4 Unsupervised CRF [6] 

1.2 MRFs and CRFs for Our Vision Problems [8] 

1.3 Our CRF Formulations for Scene Geometry Recovery Problems [10] 

1.3.1 Unsupervised CRF Formulation for Stereo Vision [10] 

1 .3.2 Unsupervised CRF Formulation for Stereo Vision Models with Monoc- 
ular Cues [12 

1.3.3 Unsupervised CRF Formulation for Depth Estimation from a Sin- 
gle Image [13] 

1 .3.4 CRF Formulation for Structure and Motion Estimation from Stereo 
Videos 

1 .4 Training Methods for Unsupervised CRFs 

1.4.1 Hard Conditional EM 

1.4.2 Contrastive Divergence [T8] 

1.4.3 Related work 

1.5 Inference using Particle-based Belief Propagation 

1.6 Applications to Scene Geometry Estimation 



VI 



Vll 

2 INFERENCE USING PARTICLE-BASED BELIEF PROPAGATION AND ITS 
APPLICATIONS FOR VISION PROBLEMS 

2.1 Belief Propagation 

2.1.1 Definitions and Notations 

2.1.2 Sum-Product versus Max-Product Algorithm 

2.1.3 Loopy BP and Message Normalization 

2.1.4 Efficiency of Updates 

2.2 Markov Chain Monte Carlo (MCMC) Methods 

2.2.1 Markov Process and Markov Chains 

2.2.2 Monte Carlo Integration 

2.2.3 Metropolis-Hasting Algorithm 

2.2.4 Simulated Annealing 

2.2.5 Gibbs sampling 

2.3 Particle Filters 

2.4 Particle-based Belief Propagation 

2.5 Applications of PBP to vision problems 

2.5.1 Particle Belief Propagation for Structure from Motion 

2.5.2 PBP for Dense Stereo Vision with Uncalibrated Cameras 0T] 

2.6 PBP for Slanted-Plane Stereo Estimation 



3 SLANTED PLANE MODEL WITH SURFACE ORIENTATION CUES FOR STEREO 
DEPTH ESTIMATION 

3.1 The Slanted Plane Stereo Model 

3.2 HOG as Surface Orientation Cues 

3.3 Relationship between HOG and Surface Orientation 

3.4 HOG Computation and Representation 

4 UNSUPERVISED CRF LEARNING FOR STEREO VISION WITH MONOCU- 
LAR CUES 

4.1 Particle Belief Propagation for Plane Estimation 

4.2 Background Review on Contrastive Divergence 

4.2.1 Introduction to Contrastive Divergence 

4.2.2 Formulation of Contrastive Divergence 

4.3 Our Unsupervised CRF Learning Algorithm 

4.4 Experimental Results 

4.5 Experimental Results with the Middlebury Dataset [77] 

4.6 Experimental Results with the Stanford Dataset [79] 

4.7 Conclusion 



Vlll 



UNSUPERVISED CRF LEARNING FOR MONOCULAR DEPTH ESTIMATION 

5.1 The model 

5.2 Stereo Pair View Prediction and Conditional Likelihood 

5.3 Unsupervised CRF Learning for MDE Model 

5.4 Experimental Results 

DEPTH AND MOTION ESTIMATION FROM STEREO SEQUENCES 

6.1 Introduction 

6.2 Proposed approach for Structure and Motion 1971 

6.2.1 Connection between Disparity and Motion 

6.2.2 3-frame Model for Depth and Motion estimation 

6.3 Experimental Results 11011 

6.3.1 Evaluation Metrics 1 101 1 

6.3.2 Results on road-driving stereo sequences 11021 

6.4 Unsupervised Learning of the Three-frame Model 11031 

6.4. 1 Unsupervised training using Hard Conditional EM/Bootstrap trainin g! 1041 

6.4.2 Unsupervised training using fourth view prediction 11071 

6.5 Conclusion fTDTI 

CONCLUSION El] 



LIST OF FIGURES 

1 . 1 Different forms of input image data for recovering scene geometry 

2.1 Particle filtering approximates the probability distribution of interest us- 
ing a large set of random samples, named particles 

2.2 The representation of the Structure from Motion problem as a bipartite 
graphical model 

2.3 Comparing bundle adjustment to PBP in structure from motion, (a) Es- 
timated mean and standard deviation of reconstruction errors for camera 
pose and map positions; (b) example posterior for one camera pose and 

(c) for one map point using PBP. |4l] 

2.4 The representation of the Dense Stereo Vision problem as a graphical 
model 

2.5 (a) Ground truth; (b) estimated depth map after a few iterations of PBP. . . 

2.6 Estimated posterior distribution of the camera pose over time 

3.1 The scene geometry is much more accurately captured by the slanted- 
plane model than by the pixel-based stereo model 

3.2 The assigned disparity plane defines a disparity value for any pixel in the 
superpixel 

3.3 HOG features for an example image 

3.4 As a surface is tilted away from the camera the edges in the direction of 
the tilt become foreshortened while the edges orthogonal to the tilt are not. 

3.5 HOG features for image regions. The amount of edge as a function of 
angle, a HOG feature, is averaged over different vertical and horizontal 
regions on various images. The different surface orientations of these 
regions affect the HOG features. We can see the cylindrical structure of 
tree trunks and the fact that the ground plane becomes more tilted as the 
distance increases 

3.6 We compute HOG features at each image scale for three different scales 
to obtain a 24-dimensional feature vector at each pixel 

4.1 Disparity plane inference algorithm using Particle-Belief Propagation. . . 

4.2 Improvement with training on the Middlebury dataset |78] 

4.3 Several depth map results on the rectified Stanford stereo dataset. The 
first 2 row are some example images. The second 2 rows are our output 
depth maps 



IX 



5.1 Resulting disparity maps from (a) the initial stereo algorithm, (b) the 
trained monocular predictor, (c) the combined monocular and binocular 
predictor 

5.2 MDE results on some testing images 



6.1 As the camera moves forward, each pixel moves along its corresponding 
epipolar line. The distance it moves depends on both their depth and 
velocity. Since the truck on the left has highest velocity w.r.t the camera, 
Pi moves the longest distance. On the other hand, the truck on the right 
has almost zero velocity w.r.t the camera, hence P<i almost stands still. 

P3 reflects the motion of the camera w.r.t the road (the scene background). 11091 

6.2 The geometric interpretation of equation (|6.1|) II 101 



6.3 . .rnn 

6.4 Top: original fourth frame In . Bottom: predicted fourth frame Id" . 
The black pixels are missing values caused by bilinear interpolation. These 
pixels are excluded when computing the view prediction error 11121 

6.5 Results on a Daimler road driving sequence 11131 



LIST OF TABLES 

4.1 Performance on the Middlebury stereo evaluation. The numbers shown 

are for unsupervised training with texture features |77] 

4.2 RMS disparity error (in pixels) and average error (average base 10 log- 
arithm of the multiplicative error) on the Stanford stereo pairs for four 
versions of our systems plus the best reported result from ( Ashutosh Sax- 



ena and Ng 2007 a[ ) on this data. Each system was either trained using 



the ground truth depth map (supervised) or trained purely from unlabeled 
stereo pairs (unsupervised) and either used texture cues (Texture) for sur- 
face orientation or did not (No texture). Note that texture information 
helps improve the performance in both supervised and unsupervised cases. 

5 . 1 Quantitative comparison between the ground plane baseline and our model 
using the average testing error per pixel as functions of the number of iter- 
ations. For each model we report: (a) The average view prediction error, 
measuring the predicted right frame and the ground truth right frame, i.e. 
the quantity on the left side of ( 5.6| ). (b) The RMS disparity error com- 



pared with the ground truth (in pixels), measuring the average difference 
per pixel in disparity value between the predicted disparity and the ground 
truth disparity 

6.1 Average view prediction error (in pixel intensity value) on 7 Daimler 

road-driving stereo sequences after each iteration of the algorithm 11031 



XI 



CHAPTER 1 
INTRODUCTION 

One of the most important characteristic of natural scenes is their underlying 3D geometric 
structure. Recovering the 3D structure of the scene from images yields a great amount 
of useful information for other computer vision tasks such as shape and scene recognition, 
object detection, or motion planning and object grasping in robotics. In this thesis, we focus 
mainly on the problem of scene geometry recovery and motion estimation. Concretely we 
want to solve the problems of monocular depth estimation, stereo vision with monocular 
cues, and structure and motion estimation from stereo video. We attempt to solve these 
problems at different levels, using different forms of input image data and different types 
of image features. 

Markov random fields (MRFs) are a very powerful tool used for modeling a wide range 
of Vision problems. Throughout this thesis, we will demonstrate that all the problems we 
are interested in can be modeled using MRFs or conditional random fields (CRFs), a variant 
of MRFs. A major part of this thesis consists of formulating general learning and inference 
algorithms for MRFs. The probabilistic inference and learning techniques we develop to 
solve our specific Vision problems can also be generalized to other problems modeled by 
MRFs and CRFs. 



1.1 Markov Random Fields and Related Formalisms 

1.1.1 MRF 

An MRF is considered a particular type of graphical model that allows us to efficiently 
compute marginal distributions. One representation of an MRF is a factor graph, or hy- 
pergraph, in which each hyperedge can connect more than one (regular) node. MRFs 
have been widely used for decades as statistical models in a variety of application areas 
Pammersley and Clifford][T97T| |Kindermann and Snell|[T980l |LifT995| ). An MRF defines 



a probability distribution on the joint assignments to (configurations of) a set of random 
variables. 

Let Y = (yi, ... , y^) be a set of iV random variables. Let Y\ , . . . , Ym be M subsets 
of Y. Let F = (/i, . . . , /m) be factors where f a is a function on assignments of values to 
the variables in Y a . 



1 M 



a=l 



where 



M 

y a=l 

Here y is an assignment of values to all random variables in Y, each y a is the projection 
of y onto the variables in the subset Y a C (y±, . . . , y^). (3 is the mode parameter. Each 
function f a is called a factor or sometimes called a local function, and is defined only on 
the subset Y a . The set of variables Y and the set of factors F define the MRF. 



3 
For inference in MRFs, the goal is to find the assignment of values to all variables 
maximizing the joint probability: 



y* = aigmax.(pp(y)) (1.3) 

V 

For training, given the i.i.d sampled training data: y , . . . , y , the objective is to maxi- 
mize the joint probability over the whole training corpus: 



T 

/3* = argmaxJJ(p j8 (2/*)) (1.4) 



P 



t=l 



The most popular version of MRFs in computer vision is the pairwise MRF model, in 
which each factor is defined over at most two random variables. 

1.1.2 CRF 

One notable variant of an MRF is a conditional random field (CRF), introduced by Laf- 
ferty et al. in (Laffe rty et al.||200T ). While MRFs model joint probabilities, CRFs model 



conditional probabilities. Similar to an MRF, a CRF is an undirected graphical model rep- 
resenting random variables and dependencies between random variables. In a CRF the 
variables are divided into two groups — exogenous and dependent. A CRF defines a con- 
ditional probability of the dependent variables given the exogenous variables and (impor- 
tantly) does not model the distribution of the exogenous variables. The advantage is that 
CRFs can estimate the distribution of interest more accurately than MRFs, since there is no 
danger of corrupting the model by modeling the exogenous variables poorly. For instance, 



4 
in the case of stereo vision one might expect that it is easier to model the probability distri- 
bution of the right image given the left image than to model a probability distribution over 
images. 

Let x = (x\, . . . xtv) represent a set of exogenous variables. In a CRF, we want to 
compute the conditional probability distribution of the dependent variables y with respect 
to the exogenous variables x. 



1 M 



where the partition function Za(x) is defined in a similar way to (1.2). 
The CRF inference equation: 



y* = &rgm&xpp(y\x) (1.6) 

y 

The CRF training equation: 

T 

P* = argmaxTT(p«(yV)) ( L7 ) 

P tl 

1.1.3 Hidden CRF 

A variation of CRFs is the hidden CRFs (HCRF) proposed by Quattoni et al. in ( |Quattoni 



et al.|2007). In HCRFs, additional latent variables Z — (z\, . . . zg) are introduced, which 



are not observed in the training data, but are part of the model. 



1 M 

Ppiv, z\x) = -TTT^ IT hfifaci) Xa, z<x) (1-8) 



where the partition function Zg(x) is defined in a similar way to (1.2). 
For inference, HCRFs predict the dependent variable y by computing p(y\x) by marginal- 
izing the hiden variable z. The HCRF inference equation: 



y* = aigmaxpp(y\x) (1.9) 

y 

= argmax^ ^pp(y,z\x) (1-10) 

V 

In training, given the training data: (x , y ), . . . , (x , y ) where x the exogenous vari- 
able, and y is the dependent variable, the objective of HCRF is to maximize the following 
conditional probability: 



T 

(3* = argmaxJTp^d/lx*) (1.11) 



P 



t=\ 



where: 



Pp(v\x) = '^2p(y,z\x) (1.12) 



1.1.4 Unsupervised CRF 

We introduce another variant of the CRF, which we call unsupervised CRF (UCRF). In 
unsupervised learning one usually formulates a parameterized probability model and seeks 
parameter values maximizing the likelihood of the unlabeled training data. Our approach 
to unsupervised learning is based on maximizing the conditional likelihood. 

The defining equation of UCRF is similar to HCRF up to some changes in variable 
naming. Here x and u denote the exogenous variables and y is the dependent variable. In 
the stereo vision example, x and u correspond to the left and right images of the stereo pair, 
while y is the hidden- state variable corresponding to the depth map of the left image. 

1 M 
P(3(u,y\x) = ^~T\ ]1 //3,a( M a,?/a,£a) (1-13) 



Zp(x 



a=l 



where the partition function Zn(x) is defined in a similar way to ( 1.2) 



For inference, UCRF aims at predicting y by computing the following: 



y* = argmaxpp(y\x 1 u) (1-14) 

y 



or 



y* = argmaxp^(?/|x) (1-15) 

y 



Note that (1.14) looks at both x and u, while ( 1.15 ) only looks at x, which in this case 
u can be considered latent variable. Later we will show that ( |1.14 ) is applied to our stereo 



inference algorithm, while ( |1.15| ) is applied to inference in our monocular depth estimator. 
In training, given the training data: (x , u l ), . . . (x , u ), the objective is to find the set 
of parameters maximizing the conditional probability of arbitrary random variables, both 
are observed by the training data. 



T 
(3* = argmax ] [ p^u^x 1 ) (1.16) 







t=l 



where: 



Pp(u\x) = ^2p(u,y\x) (1.17) 

y 

Note that in HCRF, the test time dependent variable y is part of the training data, while 
the latent variable z is not. In UCRF, the latent variable u is part of the training data, 
however the test time dependent variable y, which needs to be predicted during inference, 
is not included in the training data. 

We will then introduce a hard EM learning algorithm for maximizing the conditional 
likelihood, where we use our Particle-based Belief Propagation inference algorithm to per- 
form the hard E step, and contrastive divergence to update model parameters in the hard M 



step. Details are covered in Section [L4 



Throughout the thesis, we emphasize the employment of unsupervised CRFs learning 
for our vision models. Section [L4l will succeed this discussion with more details about the 
components of our learning approach. 

A main portion of the thesis involves applying our unsupervised CRF model for the 
scene geometry estimation problem, in particular the stereo vision problem. Here stereo 



8 
vision provides a simple setting to investigate unsupervised learning and hence seems a 
good place to start. However we believe that our approach to unsupervised learning based 
on maximizing conditional likelihood can be generalized to other, maybe much more so- 
phisticated models. 

1.2 MRFs and CRFs for Our Vision Problems 

Many problems in computer vision can be interpreted as pixel-labeling problems: the ob- 
jective is to assign labels to pixels or groups of pixels. For instance, in two-view stereo 
vision the labels can be disparity values, in segmentation the labels are indices of different 
groups that the pixels belong to, etc. 

Pixel-labeling problems are usually solved by minimizing an energy function having 
two terms: a data term penalizing assignments that are inconsistent with the observed data, 



and a smoothness term enforcing local smoothness (spatial coherence). (Geman and Ge 
man||1984 ) have pointed out that optimizing an energy function of this form is equivalent 
to estimating the maximum a posteriori probability of an MRF. 

In this thesis we mainly focus on the problem of recovery of 3D scene structure. Scene 
structure can be recovered from different forms of input image data, ranging from a single 
image to a multiple view video sequence. The list of most commonly used forms of data is 
shown in Figure [TTT] All of these image data forms are investigated in this thesis, with main 
focus on single images, stereo image pairs and stereo image pair sequences. Since images 
are 2-D representations created by perspective projection of 3-D scenes, there is a natural 
intrinsic ambiguity for recovering 3-D information from local image features. Increasing 
the amount of input image data can reduce the difficulty of the problem. 

In vision systems, depending on the form of image data used, different types of energy 
terms may arise. For the vision systems addressed in this thesis, the following types of 




Stereo Image Pair Sequence 




Stereo 
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Pair 




Monocular j: 



Image 
Sequence 




Monocular Image 









Figure 1.1: Different forms of input image data for recovering scene geometry, 
energy terms are exploited: 

• Monocular terms: depending on a single image, (e.g: intensity, gradient, SIFT, HOG, 
etc.) 

• Stereo terms: involving a stereo pair of images. 

• Motion terms: involving two consecutive video frames. 



Smoothness terms: involving only the output pixel-labeling variables. 
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As we mentioned earlier, many computer vision problems can be interpreted as pixel- 
labeling problems. The types of the output labels (dependent variables) are defined by what 
the vision problem is trying to solve specifically, i.e. the quantities of interest that we want 
to infer about the underlying scene. An image denoising algorithm would want to assign 
an intensity value or a color to each pixel, while a stereo algorithm would want to assign a 
disparity value. In this thesis, our output label types of interest are: 

• 3D point locations: for sparse Structure from Motion. 

• Disparity values: for monocular depth estimation and dense stereo depth estimation. 

• Disparity plane parameter vectors: for slanted-plane stereo model with monocular 
surface orientation cues. 

• Disparity plane parameter vectors and motion vectors: for structure and motion esti- 
mation from video sequences. 

Any vision problem we address in this thesis can be formulated by a combination of the 
subsets of each of the three elements above (input image data, types of energy terms, types 
of output labels). However, in any case, the problem can still be modeled by an MRF. 

1.3 Our CRF Formulations for Scene Geometry Recovery Problems 

1.3.1 Unsupervised CRF Formulation for Stereo Vision 

A major part of the thesis focuses on the problem of recovering scene structure from stereo 
image pairs, with emphasis on monocular surface orientation cues - evidence of surface 
orientation based on a single image. For the current discussion, let's assume we are not 
using the monocular cues. 
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A simple version of the scene structure recovery problem is dense two-view stereo 
vision, in which the depth of each pixel is recovered from a stereo image pair taken from a 
stereo camera. The depth of each pixel can be related to the horizontal disparities between 
the two stereo images. A dense output is returned - a disparity is assigned to each pixel 
in the reference image. In this problem, the exogenous variable x is the left image, the 
dependent variable y is the disparity value assignment to all pixels in x, and the exogenous 
variable u is the right image. A thorough, although not very recent, survey of state of 
the art dense stereo vision approaches was described in ( Scharstein and Szeliski|2002 ). An 



online evaluation of the best two-image stereo matching algorithms, along with test images, 
ground truth, and software can be found in ( Scharstein et al.|2001| ). 



The shift from the dense stereo model to slanted-plane stereo model started with the 
work by BirchField and Tomasi in (Birchfield a nd Tomasi| [T999). Since then, quite a few 



other stereo algorithms have been applying this slanted-plane model, among which are 
state-of-the-art algorithms in the current Middlebury benchmark (Kl aus et"aL|2006 , Wang 



and Z heng 2 008 \ |Q. Yang and Nistr||2008| ). The model uses the assumption that the scene 



is composed of a set of planar surfaces. The algorithms using this model start by segment- 
ing the image into homogeneous regions called superpixels, using image segmentation. 
The superpixels found by the segmentation method are coherent groupings of pixels hav- 
ing similar properties (color, texture, etc). Therefore it is just natural to assume that each 
superpixel corresponds to a planar surface in the scene, which from now on we will call 
disparity plane. Next, the goal of the stereo algorithm is to infer the location and orientation 
(or surface normal) of each of these planes - which are fully defined by the 3 plane parame- 
ters. After that, the disparity of each pixel can be directly computed from the disparity map 
equation, with sub-pixel precision. 

Our stereo vision model is a slanted plane model in which we use a fixed over- segmentation 
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of the left image into superpixels. The goal of inference is to assign a disparity plane for 
each superpixel. In other plane-based stereo algorithms, the authors mostly use a local 
matching step to extract reliable correspondences, followed by robust plane fitting to fit a 
disparity plane into each superpixel separately. The process is then repeated for iterative 
plane refinement. Our stereo inference algorithm also infers a disparity plane for each su- 
perpixel. However we use the Particle BP inference algorithm to simultaneously find the 
optimal joint assignment of all plane parameters to all superpixels. The optimal assignment 
is found by minimizing a high-dimensional global energy function. In this Particle BP in- 
ference algorithm, each particle, which represents a proposal disparity plane for a specific 
superpixel, gets updated over time, instead of being fixed after the plane fitting step, as in 
other slanted-plane stereo algorithms. Chapter [2] will discuss this algorithm in details. 

We formulate the problem of inferring plane parameters as an UCRF. The UCRF is a 
graphical model in which superpixels define nodes and the adjacency between superpixels 
define edges. In this case, the input image data are stereo image pairs; the energy terms 
involved are stereo and smoothness terms, and the hidden labels here are the parameters of 
disparity planes corresponding to superpixels. 

In this problem, the exogenous variable x is the left image, the dependent variable y is 
the disparity plane assignment to all superpixels in x, and the exogenous variable u is the 
right image, x and u are in the training data, while y is not. 

1.3.2 Unsupervised CRF Formulation for Stereo Vision Models with 

Monocular Cues 

Our stereo vision model is a slanted plane model involving monocular shape from texture 
cues. The introduction of an additional monocular term into the slanted-plane stereo model 
is a novel idea that has not been investigated. 
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One interesting point is that monocular cues can easily be incorporated into a stereo 
model and can be used as an additional monocular term. The monocular depth cues pre- 
dicting pixel disparity values directly can be added into a standard dense stereo model, 
without breaking the UCRF formulation of the model. For this UCRF, x and u form the 
stereo pair, and y is the assignment of disparity values to all pixels in x. Analogously, the 
surface orientation cues predicting the surface normal of superpixels can be integrated into 
our existing slanted-plane stereo model as an additional energy term. For this UCRF, x and 
u is the stereo pair, y is the assignment of disparity planes to all superpixels in x. 

Our UCRF models now involve three terms: a data matching energy (stereo term) 
measuring the degree to which the left and right images agree under the induced corre- 
spondence, a smoothness energy (smoothness term) measuring the local smoothness of the 
assigned scene structure, and a texture energy (monocular term) measuring the degree to 
which the assigned labels agree with the labels predicted by the monocular predictor. 

1.3.3 Unsupervised CRF Formulation for Depth Estimation from a 

Single Image 

Recently, there has been more interest in the problem of estimating 3-D structure from 
one single image. While this is a very interesting problem, it remains extremely more 
challenging than the stereo vision problem, due to the high ambiguity between local image 
features and the 3-D point locations, due to perspective projection. Some methods such 
as Shape from Shading ( |R. Zhang and ShahT 1999| ) make specific assumptions about the 



scene and the formation of the image, and rely mostly on photometric cues such as lighting 
and shading to recover shapes of surfaces in the scene. Several more efforts by Hoiem 
( |Derek Hoiem||2005j ) and Saxena ( |Ashutosh Saxena and Ng||2007b[ ) have demonstrated 



quite encouraging results. Both of these approaches focus on supervised learning to capture 
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the relationship between the local image features and geometric properties (orientation, 
location) of regions of the scene, as well as the relationship between different regions. In 
this thesis we use UCRF learning to train a monocular depth estimator. For MDE, we 
model the disparity of a pixel as a linear function of the monocular feature vector and a 
parameter vector. The objective of the learning algorithm is to train this monocular feature 
vector. A monocular energy term is introduced, penalizing the difference between the 
assigned disparity and the monocularly predicted disparity. An energy function formulated 
by combining this monocular energy term with a smoothness term itself defines an UCRF. 
In this UCRF, the exogenous variable x is the left image, or its feature-based representation, 
the dependent variable y is the disparity map for all pixels in x, and there is no latent 
variable u. 



For inference, we use equation (1.14), i.e, we want to infer the depth map y only by 
looking at the left image x. Training is formulated by incorporating the monocular energy 
term into a standard pixel-based UCRF stereo model. At each training step, we do infer- 
ence on this joint model to ontain the new depth map, then we use this updated depth map 
to retrain the MDE. We then implement our UCRF learning algorithm for monocular depth 
estimation for learning the monocular parameter vector using stereo image pairs only. The 
experimental results show that the monocular model is able to provide good depth predic- 
tions from different scenes in general, although there were still image parts which were not 
generalized well, i.e, it tends to perform more poorly in predicting depth for image patches 
that were too different from those in the training data. 

We also discuss in this thesis the possibility of using histograms of oriented gradients 
(HOG) feature as monocular surface orientation cues. We derive a formal relationship 
between a variant of HOG features and surface orientation. Although our observation that 
there should be a statistical relationship between HOG features and surface orientation is 
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a simple result in the area of shape from texture, HOG features have only recently gained 
popularity and to our knowledge the idea of using HOG as a surface orientation cue has not 
been previously noted. 

1.3.4 CRF Formulation for Structure and Motion Estimation from Stereo 

Videos 

Multi-view video sequence is arguably the richest form of image data for scene reconstruc- 
tion. A subtype in this form that we consider very interesting are stereo sequences. This 
data form is interesting for two reasons. First, as opposed to a multi-view video sequence, a 
stereo sequence can easily be captured using only one moving calibrated binocular camera. 
Second, this data form provides us with enough constraints to conveniently compute depth 
and motion of scene points simultaneously, since coupling dense stereo matching with mo- 
tion estimation helps decrease the number of unknowns per image point. More specifically, 
for stereo sequences, we can formulate structure and motion estimation as an energy mini- 
mization problem, in which the model is either an extension of a stereo vision model to also 
handle scene motions, or an extension of a Structure from Motion or optical flow model 
under a stereo setup. So far, to our knowledge, this useful data form has hardly been in- 
vestigated and exploited for the purpose of recovering scene structure and motion. One of 



the very few research work that has discussed this issue is (Huguet and Devernay 2007), 
in which the authors presented a variational method for scene flow estimation from a cali- 
brated stereo image sequence. In Chapter[6j we extend the constructed slanted-plane stereo 
model to handle structure and motion estimation on road-driving stereo sequences. Based 
on specific assumptions about the motion of the camera and the scene, we can reduced the 
2D optical flow problem to a ID velocity value problem. Our algorithm iteratively and 
alternately solve for structure (disparity planes) and motion (velocity values). Surface es- 
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timation is done using our slanted-plane stereo algorithm. Velocity estimation is achieved 
by solving a CRF labeling problem using Loopy BR Experiments on road-driving stereo 
sequences showed encouraging results. Performance analysis was done using our novel 
evaluation metric based on the notion of view prediction error. 

For this CRF, x is the left frame of the stereo pair at time t, u is composed of the right 
frame at time t and the left frame at time t + 1, and y is the assignment of disparity planes 
and velocity values to all superpixels in the left frame at time t. The dependent variable y 
is not known in the training video sequences. 

1.4 Training Methods for Unsupervised CRFs 



As mentioned earlier in Section 1 . 1 .4 our unsupervised CRF learning approach is based on 
maximizing conditional likelihood. This section covers more details about our learning ap- 
proach. Here we introduce a hard EM algorithm for maximizing the conditional likelihood, 
where we use contrastive divergence to update model parameters in the hard M step. 

1 .4. 1 Hard Conditional EM 



Depending on the specific use of the model, (1.17) can be further factorized as follows 
when necessary: 



Pp{u\x) = Y,Pf}(uM = J2%(y\x)Pp u Wx,y) (1.18) 



J y 
y y 



The conditional probability at the right hand side of ( 1.14) can be represented as fol- 
lows. 
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(1.19) 



For soft conditional EM, the goal of the E step is to estimate the probability distribution 
over the latent variables given the training data - this distribution can be interpreted using 



equation (1.19), while the goal of the M step is to maximize the model over complete 



data. Equation (1.19) justifies for our argument that the most likely value of the dependent 
variables given the model and the observed data, which is the result of the inference step, 
is also the value maximizing the conditional likelihood. Note that the transformation from 
(soft) conditional EM to hard conditional EM is done by replacing the sum in equation 



( 1.17 ) by the max. Hard EM works with the single most likely value of y instead of the 
distribution over y. 

Hard conditional EM is defined to be the process of iterating the following updates: 



y % := argmaxP /3 (u„2/|a; i 
V 



P 



N 



argmax V'mP«(w i ,?/ i |x l ) 



(1.20) 



(1.21) 



We will call ( 1 .20 ) the hard E step and (1.21) the hard M step. We implement the hard 
E step using our Particle-based Belief Propagation inference algorithm that we introduce 



next in Section 1.5 For the M step, we use a learning technique called contrastive diver- 



gence, more details of which will follow in Sections 1.4.2 and 4.2 These two steps are the 



foundation to the hard conditional EM algorithm that we will elaborate later in Chapter [4] 
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1.4.2 Contrastive Divergence 

In order to learn the parameters of our unsupervised CRF model, we use contrastive diver- 
gence in the hard M step of our hard EM algorithm. Contrastive divergence first introduced 
by G. Hinton et al. in ( Hinton|2002 Carreira-Perpinan and Hinton|2005 ) is a general MRF 



learning algorithm capable of training large models. Other work that also described learn- 
ing on MRF using contrastive divergence is the Field of Experts system by Roth and Black 
in dRoth and Black|2005p . 



Briefly speaking, contrastive divergence is a tool that allows us to estimate the gradient 
of a energy function, even though we cannot evaluate the energy function itself. Using con- 
trastive divergence in the hard M step allows us to compute the gradient of the energy func- 
tion with respect to the model parameters /3. We can then update /3 using gradient descent. 
By running only a few (one or two) MCMC steps, starting from the current estimation of 
the latent variable, contrastive divergence reduces significantly both the computation and 
the variance compared to running standard MCMC processes. However, since this is an 
approximate method (with finite number of samples), there is actually no known guarantee 
for convergence. 

1.4.3 Related work 



The most closely related work seems to be that of Zhang and Seitz ( |Zhang and Seitz|2007 ). 



They give a method for adapting five parameters of a stereo vision model including the 
weights for the match and smoothness energies as well as robustness parameters. The five 
parameters are tuned to each individual input stereo pair, although the method could be 
used to tune a single parameter setting over a corpus of stereo pairs. The main difference 
between their work and ours is that we train highly parameterized monocular depth cues. 
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Another difference is that we formulate a general CRF-like model for unsupervised learn- 
ing based on maximizing conditional likelihood and avoid the need for the independence 
assumptions used by Zhang and Seitz by using contrastive divergence — a general method 
for optimizing loopy CRFs ( Hinton|2002[|Carreira-Perpinan and Hinton|2005| ). 



There is also related work on learning highly parameterized monocular depth cues 
by Saxena et al. in (Ashutosh Saxena and Ng 2007b|a ), as well as Hoiem et al. in 



(Derek Hoiem] [2005 } . The main difference between these methods and ours is that we 



use unsupervised learning while they use ground truth data to train their system. One 
might argue that stereo pairs constitute supervised training of monocular depth cues. A 
standard stereo depth algorithm could be used to infer a depth map for each pair which 
could then be used in a supervised learning mode to train monocular depth cues. However, 
we demonstrate that training monocular depth cues from stereo pair data improves stereo 
depth estimation. Hence the method can be legitimately viewed as unsupervised learning 
of a stereo depth. We also demonstrate later that our stereo model with texture cues, only 
by unsupervised training, outperformed the results in ( Ashutosh Saxena and Ng|2007a| ). 



Other related work includes that of Scharstein and Pal ( |Scharstein and PaT||2007| ) and 



Kong and Tao ( Kong and Tao|[2004 ). In these cases somewhat more highly parameterized 



stereo models are trained using methods developed for general CRFs. However, the training 
uses ground truth depth data rather than unlabeled stereo pairs. 

1.5 Inference using Particle-based Belief Propagation 

Inference problems arise from various scientific areas such as Al, statistical physics, com- 
puter vision, coding theory. In MRFs,the goal of inference is to find the values of the latent 
variables achieving the maximum a posteriori of the MRF. Inference is also a vital com- 
ponent of MRFs learning algorithms in general, and of our learning approach in particular. 
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By performing inference, we want to estimate the most likely value of the latent variables 
given the model and the observed data, which can be formulated by the following equation. 



y* := aigm.axpp(y\x, u) (1-22) 

V 



From equation (1.19), we can see that (1.22) is equivalent to (1.20), which is exactly 



the hard E step in our learning algorithm described in Chapter |4j In that chapter we will 
show that the inference step minimizes an energy function defined by our MRF model. 
Recent energy minimization techniques such as Loopy Belief Propagation (LBP) ( |Pearl 



1988] |J. YedTdiaa nd Weiss 2000) or graph cuts (Y. Boy kov and Zabih||200TT ) have shown 



superior performance on MRFs than were previously possible. According to the widely 
used Middlebury stereo benchmarks ( Scharstein and Szeliski|200"2 ; Scharstein et al.|200~Tj ), 



almost all the best-performing stereo methods rely on graph cuts or LBP. However, these 
methods are known to mainly perform well on systems of many variables, each of which 
has a relatively small state space (number of possible values), or particularly nice para- 
metric forms (such as jointly Gaussian distributions). For systems with continuous valued 
variables, or discrete- valued variables with very large domains, it is impossible to directly 
use a standard discrete MRF inference techniques such as Loopy BP, graph cuts or tree- 
reweighted message passing. In our slanted plane stereo model, the MRF is a graphical 
model in which segments define nodes and the adjacency between segments define edges. 
Since the inference problem is to assign a disparity plane to segments, the state space of 
each node is a continuous 3-dimensional vector space. Similarly in our MRF model for 
Structure from Motion, our graphical model consists of two types of nodes, in which each 
camera node has 6-D state space, and each map node has a 3-D state space. For such 
systems, we propose to use a generic Particle-based Belief Propagation (PBP) algorithm 
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closely related to previous work in (Koller et al. 1999[ ), which we then formulate specif- 
ically for our MRF labeling problems. The popularity of particle filtering for inference 
in Markov chain models defined over random variables with very large or continuous do- 
mains makes it natural to consider sample-based versions of belief propagation (BP) for 
more general (tree structured or loopy) graphs. The algorithm creates an initial set of val- 
ues for each node, representing samples from the posterior marginal. It then continues 
to improve the current marginal estimate by constructing a new sampling distribution and 



draw new sets of particles. The consistency of PBP has been theoretically justified in (Ihler 



and McAllester|2009 ). The authors showed that PBP approaches the true BP messages as 



the number of samples grows large. They used concentration bounds to analyze the finite- 
sample behavior, giving a convergence rate of 0(1/ y/n) where n is the number of samples, 
independent of the domain size of the variables. Our empirical results also show that the 
algorithm is consistent in finite case, i.e. it converges to the optimal values of the message 
and belief functions with finite samples. PBP also provides estimates of uncertainty in the 
form of an approximate posterior density function. Although we only described a specific 
use of this generic PBP algorithm, we believe it can be used as an approximate inference 
scheme for a wide variety of problems that can be formulated by a probabilistic graphical 
model, especially those containing many random variables with very large or continuous 
domains. 

Similar to standard belief propagation, there are also two different versions of imple- 
mentation for PBP - sum-product and max-product. We use the max-product implemen- 
tation for the hard E step of our learning algorithm. We have also implemented the sum- 
product version for some other vision experiments that we describe later in Chapter [2] 
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1.6 Applications to Scene Geometry Estimation 

The problems we focus on in this thesis are monocular depth estimation, stereo vision with 
monocular cues, scene structure and motion estimation. We design and implement our 
own algorithms solving each of these problems. Quantitative analysis is a crucial step that 
helps evaluating the performance of our algorithms. In this thesis, we demonstrate our 
experimental results for each of our following vision systems: 

• In Chapter [2] we present our Particle-based Belief Propagation inference algorithm 
and some of its applications in vision problems. 

• In Chapter|4j we demonstrate our application of our unsupervised learning algorithm 
for stereo vision with monocular cues - the stereo model that we describe in Chap- 
ter^ The experimental results on two different datasets help validating our learning 
approach, as well as testing the performance of the final trained stereo model. 

• In Chapter |5J we describe our unsupervised CRF learning approach to monocular 
depth estimation. 

• In Chapter[6j we propose a novel evaluation metric, based on the notion of view pre- 
diction error. This metric can be used to evaluate the performance of structure and 
motion estimation algorithms on stereo sequences without ground truth data. Ex- 
perimental results on road-driving stereo sequences demonstrate that our algorithm 
successfully improves the view prediction error although it was not designed to di- 
rectly optimize this quantity. 

The thesis is organized as follows. In Chapter|2]we describe in details the application of 
our Particle-based Belief Propagation (PBP) inference algorithm to several Vision problems 
that we address throughout this thesis. Chapter [3] present our slanted-plane stereo model 
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with monocular cues, and the use of Histogram of Gradients (HOG) monocular feature 
as surface orientation cues. Chapter [4] explains our unsupervised CRF learning approach 
using hard conditional EM. Experimental results of both learning and inference on two 
different stereo datasets are also demonstrated in this chapter. We present our learning 
approach to the problem of monocular depth estimation (MDE) in Chapter [5] The unified 
framework solving for depth and motion simultaneously is introduced in Chapter [6} We 
address conclusions in Chapter|7J 



CHAPTER 2 

INFERENCE USING PARTICLE-BASED BELIEF PROPAGATION 

AND ITS APPLICATIONS FOR VISION PROBLEMS 

Although the main emphasis of the thesis is on our unsupervised CRF learning algorithm, 
one key component that takes part in both the training process and the testing process of 
our learning approach is the inference algorithm. To prepare the audience with necessary 
technical details in order to have a better understanding of unsupervised CRF learning, we 
decide to introduce our inference algorithm earlier in the thesis. 

In this chapter, we first review some background technical details based on which we 
construct our PBP inference algorithm. Next we describe in details our PBP inference 



algorithm that we mentioned earlier in Section 1.5 followed by its applications to the 
following Vision problems: 

• Structure from Motion (SfM). 

• The dense stereo vision problem with unknown camera constraints. 

• The slanted-plane stereo vision. 

By experimental results we show proofs of the convergence for the accuracy of such a 
generic PBP algorithm with finite samples. In addition to accuracy, the PBP algorithm also 
allows us to estimate and represent state uncertainty, although at some computational cost. 
In particular, an experiment with a synthetic structure from motion arrangement shows that 
its accuracy is comparable with the state-of-the-art Bundle Adjustment (BA) while allowing 
estimates of state uncertainty in the form of an approximate posterior density function. 

24 



25 

2.1 Belief Propagation 

2.1.1 Definitions and Notations 

Let G be an undirected graph consisting of nodes V = {1, . . . , k} and edges E, and let 
T s denote the set of neighbors of node s in G, i.e., the set of nodes t such that {s, t} is 
an edge of G. In a probabilistic graphical model, each node s E V is associated with a 
random variable X s taking on values in some domain, X s . We assume that each node s 
and edge {s, t} are associated with potential functions ty s and \I/ S t respectively, and given 
these potential functions we define a probability distribution over assignments of values to 
nodes as 



Mn **(*-)) ( n *-,*#»» 2*)) 

V s J \{ s ,t}eE ) 



P(x) = ^\ \ *«(**) *«,t(a?*,a:*) (2-1) 



Note that (2.1 ) corresponds to the pairwise MRF model. Here x is an assignment of 
values to all k variables, x s is the value assigned to X s by x , and Z is a scalar chosen to 
normalize the distribution P (also called the partition function). We consider the problem 
of computing marginal probabilities, defined by 

p s (x s )= Yl p ^ < 2 - 2 ) 

X . X g — Xg 

In the case where G is a tree and the sets X s are small, the marginal probabilities can 



be computed efficiently by belief propagation (Pearl 1988). This is done by computing 



messages mt^s each of which is a function on the state space of the target node, X s . 
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These messages can be defined recursively as 

mt-> a (x a )= JZ ^t,s( x t' x s)^t(xt) JJ m u ^ t (x t ) (2.3) 

When G is a tree this recursion is well founded (loop-free) and the above equation uniquely 
determines the messages. We will use an unnormalized belief function defined as follows. 

B s (x s ) = V s (x s ) H mt^ s (xs) (2.4) 

ter s 

When G is a tree the belief function is proportional to the marginal probability P s defined 



by (2.2). It is sometimes useful to define the "pre-message" Mt^ s as follows for xt € Xf. 



M t ^ s {x t ) = ^t(xt) JJ m u ^ t {x t ) (2.5) 

ueT t \s 

Note that the pre-message M^ s defines a weighting on the state space of the source node 
Xt, while the message mt^ s defines a weighting on the state space of the destination, X s . 



We can then re-express (2.3 )-(2.4) as 



mt^ s (x s ) = 2_^ ^t 1 s{xt,x s )Mt^ s (x t ) Bt(xt) = Mt-> a (xt)m a -n(,xt) 

x t eX t 

Although we develop our results for tree-structured graphs, it is common to apply belief 
propagation to graphs with cycles as well ("loopy" belief propagation). We note connec- 
tions to and differences for loopy BP in the text where appropriate. 

For reasons of numerical stability, it is common to normalize each message mt^ s so 
that it has unit sum. However, such normalization of messages has no other effect on 



the (normalized) belief functions (2.4). Thus for conceptual simplicity in developing and 
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analyzing particle belief propagation we avoid any explicit normalization of the messages; 
such normalization can be included in the algorithms in practice. 

Additionally, for reasons of computational efficiency it is common to use the alternative 
expression mt-> s { x s) — S ^t s( x t, x s)Bt(xt) / m s^t( x t) when computing the messages. 
By storing and updating the belief values Bt(xt) incrementally as incoming messages are 
re-computed, one can significantly reduce the number of operations required. Although 



our development of particle belief propagation uses the update form (J23J), this alternative 
formulation can be applied to improve its efficiency as well. 

2. 1 .2 Sum-Product versus Max-Product Algorithm 

There are two different algorithms that can be chosen to implement BR The sum-product 
algorithm computes the marginal distribution at each node of the MRF. The max-product 
algorithm, on the other hand, aims at computing the MAP estimate of the whole MRF. 
Therefore, max-product belief propagation attempt to solve the same problem as that of 
Graph Cuts algorithms. 

2.1.3 Loopy BP and Message Normalization 

Loopy belief propagation maintains a state which stores a numerical value for each mes- 
sage value mt-> s (x s ). In loopy BP one repeatedly updates one message at a time. More 
specifically, one selects a directed edge t — V s and updates all values mt^. s (x s ) using 
equation ( |2.3[ ). Loopy BP often involves a large number of message updates. As the num- 



ber of updates increases message values typically diverge — usually tending toward zero 
but possibly toward infinity if the potentials ty s and ^!t,s arQ large. This typically results 
in floating point underflow or overflow. To avoid floating point errors messages can be 
normalized — the function mt^. s can be scaled by a constant so that it sums to one. We 
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can normalize an entire state by normalizing each message. It is important to note that 
normalization commutes with update. More specifically, normalizing a message, updating, 
and then normalizing the resulting state, is the same as updating (without normalizing) and 
then normalizing the resulting state. This implies that any normalized state of the system 
can be viewed as the result of running unnormalized updates and then normalizing the re- 
sulting state only at the end. For conceptual simplicity we avoid explicit normalization 
throughout this paper. But all algorithms described here can be normalized and have the 
property that message updates commute with normalization in the sense just mentioned. 

2. 1 .4 Efficiency of Updates 



Performing a message update using equation (2.3) would seem to involve an inner loop 



over the node u. This inner loop can be avoided by using the following which is equivalent 



to (2.3). 



The belief values Bt can be stored as part of the state and maintained incrementally when 
updates are made. This efficiency improvement is possible with the particle belief algo- 
rithm described in the next section. However, for conceptual simplicity we will consider 



only the inefficient update analogous to (2.3 ). 



2.2 Markov Chain Monte Carlo (MCMC) Methods 

Markov chain Monte Carlo (MCMC) is a class of algorithms for simulating direct draws 
from complex probability distributions of interest. MCMC has its name from the fact that 
it uses the previous samples to randomly generate the next samples, constructing a Markov 
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chain that has the desired distribution as its equilibrium distribution. The state of the chain 
after a large number of steps is then used as a sample from the desired distribution. The 
quality of the sample improves as a function of the number of steps. 

2.2.1 Markov Process and Markov Chains 

Let Xt denote the value of a random variable at time t, and let the state space refer to the 
range of possible X values. The random variable is a Markov process if the transition prob- 
abilities between different values in the state space depend only on the random variables 
current state, i.e., 



PriX t+1 \X t ,X t ^ 1 ,...,X ) = Pr(X t+1 \X t ) (2.7) 

Thus for a Markov random variable the only information about the past needed to pre- 
dict the future is the current state of the random variable, knowledge of the values of earlier 
states do not change the transition probability. A Markov chain refers to a sequence of 
random variables (Xq; ...; X n ) generated by a Markov process. 

2.2.2 Monte Carlo Integration 

The Monte Carlo approach was originally used as a method to compute integrals approxi- 
mately by random number generation. Suppose we wish to compute a complex integral that 
can be represented as a product of a function f(x) and a probability density function p(x), 
the integral can then be expressed as an expectation of f(x) over p(x). This expectation in 
turn can be approximated by drawing a large number of random variables fromp(x): 
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rb rb i n 

/ g(x)dx = f(x)p(x)dx = E p ( x )[f(x)] « - ^ /(xj (2.8) 

J a J a ' n 

2.2.3 Metropolis-Hasting Algorithm 

The most challenging problem of applying Monte Carlo integration is how to sample from 
some complex probability distribution p(x). Suppose our goal is to draw samples from 
some distribution p(x) where p(x) = j?f(x), where the normalizing constant Z may not 
be known, and very difficult to compute. 

The Metropolis algorithm (Metropol is and Ulam||1949 ) generates a sequence of draws 



from this distribution by first computing 



. ( f(9*)Q(9 t -i,9n \ 
\f(0t-l)Q(0*,8t-l) J 

at each time step t, and then accepting a candidate point with probability a. If the 
candidate is rejected, the current value 9f_\ is retained. Here Q(9t-i, &*) is the transition 
probability function. The original Metropolis algorithm calls for the transition probability 



function to be symmetric, i.e. Q(6t-i, 0*) = Q(6*,6t-i); Hastings in (Hastings) general 



ized the Metropolis algorithm by using an arbitrary transition probability function. 

2.2.4 Simulated Annealing 

Simulated annealing was developed as an approach for locating a good approximation to 
the global optimum of a complex multimodal function in a large search space. This is a 
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good solution for cases where standard hill-climbing approaches such as gradient descent 
may get stuck at a local optimal peak. 

The idea is that when we initially start sampling the space, we will accept a reason- 
able probability of a down-hill move in order to explore the entire space. As the process 
proceeds, we decrease the probability of such down-hill moves. 



— -''sH^M' 



We can see that simulated annealing is very closely related to Metropolis sampling, 
differing only in that the probability a of a move also depends on a temperature T(t), often 
called the cooling schedule. If T is set to 1 then equation ( 2.10[ ) is exactly ( |2.9[ ). Simulated 



annealing usually starts off with a high temperature - high transition probability, and then 
cool down to a very low temperature (T m 0) - low transition probability. 

2.2.5 Gibbs sampling 

Gibbs sampling is considered a special case of Metropolis-Hastings sampling wherein the 
random value is always accepted (i.e, a = 1). The purpose of Gibbs sampling is to generate 
a sequence of samples from the joint probability distribution of multiple random variables. 
The purpose of such a sequence is to approximate the joint distribution, or to compute an 
integral (such as an expected value). 

The key to the Gibbs sampler is that one only considers univariate conditional distribu- 
tions the distribution when all of the random variables but one are assigned fixed values. 
Such conditional distributions are far easier to simulate than complex joint distributions and 
usually have simple forms. Thus, one simulates n random variables sequentially from the 
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n univariate conditionals rather than generating a single n-dimensional vector in a single 
pass using the full joint distribution. 

2.3 Particle Filters 

Particle filtering, also known as sequential Monte Carlo methods (SMC), are sophisticated 
model estimation techniques based on simulation (( Arulampalam et al.|2002| Doucet et al 
200T||Godsill and Clapp|200T||Khan et al.|2005| )). 



They are usually used to estimate Bayesian models and are the sequential ('on-line') 
analogue of Markov chain Monte Carlo (MCMC) batch methods and are often similar to 
importance sampling methods. Well-designed particle filters can often be much faster than 
MCMC. They are often an alternative to the Extended Kalman filter (EKF) or Unscented 
Kalman filter (UKF) (( van der Merwe et al.|2001' )) with the advantage that, with sufficient 



samples, they approach the Bayesian optimal estimate, so they can be made more accurate 
than either the EKF or UKF. However, when the sample size is not sufficiently large, they 
might suffer from sample impoverishment. Particle filters can also be combined by using a 
version of the Kalman filter as a proposal distribution for the particle filter. 

Particle filtering starts with a large set of random samples called particles and their 
weights: 



{(xiwl):i = l,..,p) 



P 



5>i = i 



i=l 
Particle filtering then iteratively approximates the probability distribution of interest 

using the initial set of particles. 
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Figure 2.1: Particle filtering approximates the probability distribution of interest using a 
large set of random samples, named particles. 



f(x t )p(x t \yo, .., yt)dxt « J^ «4/(a;|) 

i=l 



A single iteration of the algorithm (as illustrated in Figure 2.1 lis as follows 



1 . Draw samples from the proposal distribution. 



2. Update the importance weights up to a normalizing constant. 



3. Draw P new particles from the current particle set with probabilities proportional to 
their weights. 
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4. Set the weight of each particle to be 1/P. Update the proposal distribution using the 
new set of particles. 

2.4 Particle-based Belief Propagation 



We now consider the case where \X S \ is too large to enumerate in practice and define a 
generic particle (sample) based algorithm. This algorithm essentially corresponds to a non- 
iterative version of the method described in ( |Koller et al.|1999[ ). 



The procedure samples a set of particles x s , . . ., x s with x s G X s at each node s of 
the networl|^] drawn from a sampling distribution (or weighting) W s (x s ) > (correspond- 
ing to the proposal distribution in particle filtering). First we note that ( |2.3[ ) can be written 
as the following importance- sampling corrected expectation. 



mt->s{x s ) 



E 



Xt ~W t 



^ S: t(x s ,x t )^t(xt) Yl U £T t \s mu^t{xt) 
W t (x t ) 



(2.11) 



Given a sample x\ , . . ., x t n of points drawn from Wt we can estimate mt-> s (x s ) as 



follows where wf' = W s (xg). 



m 



(i) 



1 - 
n ^ 



Hs(x?\^ ) )M4 J) )Tlu£r t \ s ™ { jl>t 



j=i 



w 



U) 



(2.12) 



Equation ( |2.12| ) represents a finite sample estimate for ( |2. 1 1 [ ). Alternatively, ( |2.12[ ) defines a 
belief propagation algorithm where messages are defined on particles rather than the entire 



1. It is also possible to sample a set of particles {xf t } for each message m s ^ t in the net- 
work from potentially different distributions W s t(x s ), for which our analysis remains essentially 
unchanged. However, for notational simplicity and to be able to apply the more computationally 
efficient message expression described in Section |2.1| we use a single distribution and sample set 
for each node. 
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set X s . As in classical belief propagation, for tree structured graphs and fixed particle 
locations there is a unique set of messages satisfying ( |2.12[ ). Equation ( |2.12| ) can also be 

applied for loopy graphs (again observing that message normalization can be conceptually 

(i) (i) 

ignored). In this simple version, the sample values x s and weights w s remain unchanged 

as messages are updated. 



We now show that equation ( |2.12[ ) is consistent — it agrees with ( |2.3[ ) in the limit as 
n — >■ oo. For any finite sample, define the particle domain X s and the count Cj(x) for 
x G X s as 



X s = {x s eX s : 3k x s l) =x s } 



(i) 
Cs{x s ) = \{k: x s ' =x s }\ 



Equation (2.12) has the property that if x 



(0 



x; s then m 



t-ts 



m 



ii') . 



t—ts 



; thus we can 



rewrite (2.12) as 



mt^s{Xsj 






" t -^ W t (x t ) 
x t £Xt 



Vt,s{xux a )* t (xt) n fh u ^ t {x t ) x s eX s (2.13) 
u€T t \s 



Since we have assumed W s (x s ) > 0, in the limit of an infinite sample Xf becomes all of 



Xt and the ratio (ct(xt)/n) converges to Wt(xt). So for sufficiently large samples ( |2.13| ) 
approaches the true message ( |2.3| ). Fundamentally, we are interested in particle-based ap- 
proximations to belief propagation for their finite-sample behavior, i.e., we hope that a rel- 
atively small collection of samples will give rise to an accurate estimate of the beliefs - the 
true marginal Pf(x-t). At any stage of BP we can use our current marginal estimate to con- 

(i) 

struct a new sampling distribution for node t and draw a new set of particles {x\ }. This 
leads to an iterative algorithm which continues to improve its estimates as the sampling 
distributions become more accurately targeted. Unfortunately, such iterative resampling 
processes often require more work to analyze; see e.g. ( Neal et al.|2003 ). 
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Map points 



Camera poses 




Figure 2.2: The representation of the Structure from Motion problem as a bipartite graphi- 
cal model. 



In (Koller et al. 1999), the sampling distributions were constructed using a density 
estimation step (fitting mixtures of Gaussians). However, the fact that the belief estimate 
Mt(xt) can be computed at any value of xt allows us to use another approach, which has 
also been applied to particle filters ( |Godsill and Clapp|2001||Khan et al.|2005| ) with success. 
By running a short MCMC simulation such as the Metropolis-Hastings algorithm, one 
can attempt to draw samples directly from the weighting Mp This manages to avoid any 
distributional assumptions inherent in density estimation methods, but has the disadvantage 
that it can be difficult to assess convergence during MCMC sampling. 



2.5 Applications of PBP to vision problems 

In this section we describe the applications of the PBP algorithm to several Vision problems 
related to scene geometry recovery. 



2.5.1 Particle Belief Propagation for Structure from Motion 

In structure from motion(SfM), given a set of 2D images of the same scene, the objective 
is to simultaneously recover both the camera trajectories and the 3D structure of the scene. 
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For this problem, it is commonly accepted that Bundle Adjustment (BA) (|Triggs et al.|1999|) 



is the Gold standard. There has been related work on this problem (Hartley and Zisserman 



2004||Kaucic et al.|2001||Snavery et al.|2006[ ), mostly using geometric based methods. All 



of these works demonstrated that their result using BA as the last step is much better than 
without BA. With each camera pose and each map object 3D location being represented 
as a finite set of particles, we show that PBP can successfully estimate their true states by 
estimating their posterior distributions over the state space given the image observations. 
The algorithm was tested on both real and synthetic data. An experiment on synthetic data 
allows us to compare the performance of our method with BA. 

First, we represent the problem of SfM in the context of PBP. Each observed image is 
converted to a sparse set of image points. These points can be high-level image features, 



obtained by a feature detector (such as corner detector or SIFT detector (Lowe 2004)). 
The correspondences of image points between images are automatically obtained using a 
feature matching method. This also defines the correspondences between the image points 
and the map points. The scene is represented as a sparse map of 3D points. Each camera 
pose is also a 3D point, combined with 3 angles of rotation, which define the rotation of 
the camera about 3 axes. Our method is based on the assumption that given a set of image 
observations and the estimate of all map points, there exists a probability distribution for 
each camera pose over its state space. The true state of map points and camera poses are 
hidden variables that we want to estimate. Let Pj e 3ft denote the random variable for 
the i camera pose. Let Yj E 3ft denote the random variable for the j map point. The 
observed data are image points and their correspondences. We denote x^j G 3ft the image 
point variable associated with map point Yj as seen from camera P^. 

Our graphical model G(V,E) consists of two types of nodes. Each camera node % is 
associated with a camera variable Pj, each map node j is associated with a map variable 
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Yj. For simplicity of notation we name each node after its variable. For each camera and 
a map point it observes, we add an edge connecting two respective nodes into the graph. 
No edge connects any two camera nodes, or any two map nodes. If each map point is seen 
by all cameras, we have a complete bipartite graph. For each edge in the graph, there is a 
binary potential function associated to it, denoted \E^ j(P^, Yj). More specifically, we have: 



II l|2 

-\\0(P- Y-)— x- • 

VijiPi, Yj) = e 2*^ (2.14) 



where Q(P{, Yj) is the reprojection function that takes a camera pose and a map point 
and returns the reprojected image point, x^j is the image observation of map point yj as 
seen from camera pi, o is the variance of the Gaussian image noise. This is equivalent to 
assuming a normal distribution of the true observation around the given observation. This 
allows us to represent the uncertainty of the observations. 



At this point the message function and belief function on G are well defined using ( 2.3 ) 



and ( |2.4[ ). However, the state space of each hidden variable in G is too large to do inference 
with BR In order to use PBP, we discretize the state space into a relatively small number 
of states, each of which is represented by a particle. These states can be sampled from 
a normal distribution with some initial mean and variance. Now we have at each camera 
node Pi, a set of M particles: P. , . . ., P^ , at each map node Yj, a set of N particles: 
jn ' , . . ., y- . If we assume no unary potential at each node, we can define the message 



going from camera node i to map node j in G by (2.12) as follows: 



i M 

^-iS^^^W") II 41: (2-15) 

h=i uei\\j 



m 
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At the beginning, all particles in a node are assigned uniform weights and no message 



has been computed, equation (2.15) becomes: 



M 

«2li = s^*« ( ^W } <2 ' 16) 

h=l 

This can be interpreted as the marginal probability J2h=l P ( Y j = ^ |sy,Pi = 
P- ). It follows that the belief of a map node becomes the marginal probability over 
all assignments of its neighboring camera nodes, conditioned on relevant observations: 



M 

B J (Yj k) )=m2 p ( Y J= Y J k) \ x > p *= p t ) ) =Y, p t Y j= Y J k) \ x ' p ) < 2 - 17 ) 

i€Tj h=l p 



where P is an assignment of values to all neighboring pose nodes. The message from 
a map node to a pose node and the belief of a pose node can be interpreted similarly. As 
BP proceeds, the information from one node is sent to all other nodes in the graph. This 



allows nodes of the same type to communicate with each other. Then (2.17) will become 



exactly (2.2), which is what we want to compute. As shown in (Ihler and McAllester 



2009), the posterior marginal estimated by PBP is in most cases guaranteed to converge to 
the true posterior marginal distribution estimated by BP. In addition to PBP, the samples at 
each node are iteratively updated by Gibbs sampling from the estimated posterior marginal 
distribution at that node. Gibbs sampling allows particles to freely explore the state space 
and thus compensates the inadequacy of representing a large state space by a small number 
of samples. After a sufficient number of alternative PBP and Gibbs sampling iterations, the 
estimated posterior at each node converges to the true posterior distribution, which gives us 
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the answer for the final state of each camera pose and map point. 

The objective of this experiment is to compare PBP with the well known method BA 
in terms of reconstruction accuracy, by measuring the deviation of their reconstruction 
estimates from the ground truth. The model for bundle adjustment is a pairwise graphical 
model over poses pi and object positions yj given observed locations X{j in the images: 



II II 2 

-||Q(Pi>2/j)-Zjj|| 

p ({Pi}, {yj}) = Yi^ijiPhVj) ^i,j(Pi,Vj) = e 2^2 



Bundle adjustment then searches for a local maximum (i.e., a mode) of the posterior 
distribution over the {pi, yj}. However, if the intent is to minimize the expected squared 
error loss, we should prefer to estimate the mean of the posterior distribution rather than 
its mode. These can be very different in problems where the posterior distribution is very 
skewed or multi-modal; in such cases, it may be advantageous to estimate the posterior dis- 
tribution and its mean using methods such as belief propagation. Additionally, an explicit 
representation of uncertainty can be used to assess confidence in the estimate. 

We show a comparison between estimates given by BAnand PBP on synthetic structure 



from motion data in Figure 2.3 For the data, we generate a fixed set of map points and a 
few camera poses. Each point's observations are obtained by projecting the point on each 
camera's image planes and adding Gaussian noise. We assume that the camera calibration 
parameters and the image correspondences are known, and initialize the estimates for both 
algorithms to ground truth, then run each to convergence and compare the resulting error. 



Figure 2.3 a) shows the estimated mean and standard deviation of reconstruction errors 



of BA and PBP from 200 different runs (units are synthetic coordinates). This shows that, 



2. We use the sba package in (Lourakis and Argyros 



2004). 
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Method 


Error 


BA (Map) 


185.81 ±66.43 


PBP (Map) 


194.00 ± 63.02 


BA (Pose) 


11.88 ±4.40 


PBP (Pose) 


11.56 ±4.39 



Camera Pose Posteror Estimate 



.>.'■. 



Map Point Posteror Estimate 



(a) 



(b) 



(c) 



Figure 2.3: Comparing bundle adjustment to PBP in structure from motion, (a) Estimated 
mean and standard deviation of reconstruction errors for camera pose and map positions; 
(b) example posterior for one camera pose and (c) for one map point using PBP. 

at the same level of image noise and from the same initial conditions, PBP produces essen- 
tially the same accuracy as BA for both camera and map points. However, we expect that in 
less idealized cases (including, for example, incorrect feature associations or outliers mea- 
surements), PBP may perform much better. In addition, PBP provides an estimate of state 
uncertainty (although at some computational cost). Figure |23| b)-(c) shows the estimated 
posterior distributions given by PBP for a camera pose and map point, respectively, along 
with the mean found via PBP (circle), mode found via BA (diamond), and true position 
(square). 



2.5.2 PBP for Dense Stereo Vision with Uncalibrated Cameras 

Classical dense two-frame Stereo algorithms compute a dense disparity map for each pixel 
given a pair of images under known camera constraints (i.e, the orientation of the 2 cameras 
and the distance between them are known) ( |Sun et al.|2003[[Felzenszwalb and Huttenlocher 



2006, Scharstein et al. 2001). In this experiment, given a pair of stereo images with un- 
known camera constraints, we use PBP to simultaneously compute the dense depth map 
for each pixel, and the configuration of the second camera relative to the first. 

The formulation of the graphical model in this case is quite similar to the previous 
problem. However there are only 2 camera nodes, and the number of map nodes equals the 
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Figure 2.4: The representation of the Dense Stereo Vision problem as a graphical model. 
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Figure 2.5: (a) Ground truth; (b) estimated depth map after a few iterations of PBP. 

number of pixels in the first image. An edge is added between every pixel and every camera. 



(see Figure 2.4) Figure 4.2 shows the estimated depth map after a few iterations. We do 
not show evaluation results of our output depth map in comparison with ground truth and 



other stereo algorithm (using for example the Middlebury stereo benchmark (Scharstein 
et al.|200T )), because such a comparison will not be very meaningful, as the labels we use 
are the true 3D depths of each pixel instead of their disparities. However, in Figure [2^6] we 



plot the estimated posterior distribution of the second camera pose at each iteration, thus 
show that the distribution gradually approaches the true state over time. 
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Figure 2.6: Estimated posterior distribution of the camera pose over time 



2.6 PBP for Slanted-Plane Stereo Estimation 



For this task, given a pair of images (X, Y), and a given segmentation of X, and a given 
setting of the model parameters, the inference problem is to find an assignment Z of plane 



parameters to superpixels so as to minimize the total energy in (3.2). The energy defines a 
Markov random field. More specifically, the match energy defines the unary potential on 
each superpixel independently and the smoothness energy defines the binary potential on 
pairs of adjacent superpixels. 



The complete inference algorithm is described in Section 4.1 



CHAPTER 3 

SLANTED PLANE MODEL WITH SURFACE ORIENTATION 

CUES FOR STEREO DEPTH ESTIMATION 

In this chapter, we present our slanted-plane stereo vision model with monocular cues for 
surface orientation estimation. We also explain our idea of using HOG features as surface 
orientation cues, and present an analytical justification for this idea. 

3.1 The Slanted Plane Stereo Model 

As previously mentioned in our introduction, the two-view stereo vision problem is com- 
monly treated as a pixel-labeling problem, where labels are disparity values, usually quan- 
tized to be integers. The disparity is equivalent to the horizontal displacement of the pixel 
between the left and the right image, or the inverse depth of that pixel. Specifically, the 
corresponding pixel in the right image of a pixel p = (x p , y p ) in the left image is given 
by p' = p — d(p) = (x p — d(p), yp). The labeling is then resolved by using an energy 
minimization technique. The standard form of the global energy function is formulated in 



equation (3.1 ) 



E{d) = E M {d) + \E s {d) (3.1) 

where Ej^ is the match energy measuring how well the disparity assignment d agrees 
with the input image pair, Eg is the smoothness energy enforcing the smoothness assump- 
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Figure 3.1: The scene geometry is much more accurately captured by the slanted-plane 
model than by the pixel-based stereo model. 

tion. Stereo methods aiming at finding d that optimizes the global energy function are called 



global methods. The best-performing algorithms in the Middlebury benchmark ( Scharstein 
et al.|200T ) are all global methods. 

The quantization of disparity values d(p) helps reduce considerably the state space for 



searching. All techniques discussed in (Scharstein and Szeliski 2002) follow a common 
framework - producing the quantized disparity map as output. However, this output has the 
effect that the scene seems to be split into a series of fronto-parallel surfaces, as illustrated 



in Figure 3.1 In other words, piecewise continuity of the scene is replaced by piecewise 



constancy. This is a poor explanation of the scene geometry, especially when the scene 
contains slanted surfaces. 

The use of slanted-plane model for stereo vision was first proposed by BirchField and 
Tomasi in ( Birchfield and Tomasi||1999 ) as a special case of the affine model. Since then 
the model has been widely applied in many other stereo algorithms, including most of 
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state-of-the-art algorithms in the current Middlebury benchmark. ( Klaus et al.|2006 ; Wang 



and Zheng 2008, Q. Yang and Nistr 2008). The model assumes that the scene structure 



is locally planar - the scene is composed of a set of planar regions. This is a reasonable 
assumption, based on the fact that most surfaces found in a natural environment can be 
approximated by a plane or a set of planes. These algorithms start by segmenting the image 
into homogeneous regions called superpixels, using an image segmentation method. The 
superpixels found by the segmentation method are coherent groupings of pixels having 
similar properties (e.g: color, texture). Therefore it is just natural to assume that each 
superpixel corresponds to a disparity plane. Next, the goal of the stereo algorithm is to 
infer the location and orientation (or surface normal) of each of these planes - which are 
fully defined by the 3 plane parameters. After that, the disparity of each pixel can be 
directly computed from the disparity map equation, with sub-pixel precision. 

Our stereo vision model is a slanted plane model involving monocular shape from tex- 
ture cues. The introduction of an additional monocular term into the slanted-plane stereo 
model is a novel idea that has not been investigated. Another difference of our work is in 
the inference algorithm. Other related work that also employ the plane-based stereo model 
use different approach in the inference step. One of the best methods in the Middlebury 
benchmark is (Klaus et al. 2006| ). The authors used a local matching step to extract re- 
liable correspondences, followed by robust plane fitting to fit a disparity plane into each 
superpixel separately. The process were repeated for iterative plane refinement. For plane 
fitting, instead of trying to estimate all plane parameters jointly, the authors proposed a 
decomposition method to solve each parameter, again separately. The estimated planes are 



fixed after this step. Other slanted-plane stereo algorithms as in (Wang and Zheng 2008 



Q. Yang and Nistr|20 08 ) also follow a very similar approach for plane estimation. 



Our stereo inference algorithm also infers a disparity plane for each superpixel. How- 
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ever we use the Particle BP inference algorithm to simultaneously find the optimal joint 
assignment of all plane parameters to all superpixels. The optimal assignment is found 
by minimizing a high-dimensional global energy function. With this PBP inference algo- 
rithm, each particle, which represents a proposal disparity plane for a specific superpixel, 
gets updated over time, instead of being fixed after the plane fitting step, as in previously 
mentioned methods. 

The energy function involves three terms — a correspondence energy measuring the 
degree to which the left and right images agree under the induced disparity map, a smooth- 
ness energy measuring the smoothness of the induced depth map, and a texture energy 
measuring the degree to which the surface orientation at each point agrees with a certain 
(monocular) texture based surface orientation cue. Specifically, the global energy in equa- 



tion (3.1 ) is replaced by the following: 



E{Z) = E M {Z) + E S (Z) + E T {Z) (3.2) 

where Z is the assignment of a disparity plane to each superpixel in the left image. We 
denote X to be the left image, Y to be the right image. For each superpixel i of X we have 
that Z specifies three plane parameters A^ B^, and C{. Given an assignment Z of plane 
parameters to superpixels we define the disparity d(p) for any pixel p as follows (where 
i(jp) is the superpixel containing p and x p and y p are the image coordinates of p). 



d (p) = A i{ P ) x P + B t( P )Vp + C i(p) (3-3) 



So by equation (3.3 ) we have that Z assigns a real valued disparity to each pixel. 
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Figure 3.2: The assigned disparity plane defines a disparity value for any pixel in the su- 
perpixel. 

To define the smoothness energy we write (j>, q) E Bjj if p is a pixel in superpixel i, q 
is a pixel in superpixel j, and p and q are adjacent pixels (p is directly above, below, left or 
right of q). The smoothness energy is defined as follows (where t$ and \g are parameters 
of the energy, i ranges over all superpixels, N(i) is the set of superpixels adjacent to i, and 
z is the overall assignment of three dimensional vectors of plane parameters {A\, B^, Q) 
for all superpixels i): 



( 



E S (Z) 



y min 

tJeJV(i) 



\ 



V 



r s , Yl *s\d(p) - d(q)\ 



(3.4) 



/ 



Intuitively the minimization with rg corresponds to interpreting the entire boundary 
between i and j as either an occlusion boundary or as a joining of two planes on the same 
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object. 

Next we consider the match energy. We write p — d(p) for the pixel in Y that corre- 
sponds to the pixel p in image X under the disparity d(p). For color images we construct 
a nine dimensional feature vector $ (p) and $ (p) for the pixel p in the images X and 
Y respectively. The vector $ (p) consists of three (bias gain corrected) color values plus 
a six dimensional color gradient vector and similarly for $p . We write $^ (p) for the kth 
component of the vector $ (p). The match energy is defined as follows where A^ are nine 
scalar parameters of the match energy. 

e m(z) = J2 E x ^t (p) - ®i(p + d (p)) ) 2 ( 3 - 5 ) 

P eX k 

Finally we consider the texture energy. At each pixel p we also compute a HOG vector 
H(p) which is a 24 dimensional feature vector consisting of three 8 dimensional normal- 
ized edge orientation histograms — an 8 dimensional orientation histogram is computed 
at three different scales. The texture energy is defined as follows where i(p) is the super- 
pixel containing pixel p and where the scalars tt, A^, A#, and the vectors ft a and fig are 



parameters of the energy. The form of this energy is justified in Chapter 3.2 



^ I A A (d(p)((3 A -H(p))-A i{p) ) \ 

p \ + A B (d{p){(l B -H{p))-B i{) ) 



(3.6) 



These energy terms can also be interpreted in a probabilistic way as follows. The com- 
bination of the smoothness energy term and the texture energy term can be interpreted as an 
energy E%(X, Z, f3z), which determines the probability P(Z\X, fiz)- The match energy 
term can be interpreted as an energy Ey(X, Y, Z, /3y), which determines P(Y\X, Z, /3y), 
The general conditional probability formulation of our model will be discussed in details 
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in Chapter |4| 



3.2 HOG as Surface Orientation Cues 
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Figure 3.3: HOG features for an example image. 



3.3 Relationship between HOG and Surface Orientation 



In this chapter we discuss the possibility of using a specific type of monocular cues - the 
HOG feature (Histogram of Oriented Gradients). This feature has recently been exten- 
sively applied in computer vision and image processing for the purpose of object detection 



((Dalai and Triggs 2005b, Felzenszwalb et al. 2008)). The technique counts occurrences 



of gradient orientation in localized portions of an image. This method is similar to that of 



edge orientation histograms, scale-invariant feature transform descriptors (or SIFT (Lowe 
2004)), and shape contexts ( |S. Belongie and Puzicha||2002[ ), but differs in that it computes 
on a dense grid of uniformly spaced cells and uses overlapping local contrast normalization 
for improved performance. 
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Image Plane 




Figure 3.4: As a surface is tilted away from the camera the edges in the direction of the tilt 
become foreshortened while the edges orthogonal to the tilt are not. 

First we observe that HOG, which describe the edge distribution, can be related to the 
orientation of surfaces. We formulate the statistical relationship between HOG features and 
surface normal. Based on this derivation, we then argue that HOG features can be used as 
surface orientation cues, and thus can be incorporated in our existing slanted-plane stereo 



model (as described in Section [371] ) . 

Here we aim at describing a justification for the idea of using HOG features as surface 
orientation cues, as mentioned in Chapter[3j The basic intuition behind HOG as an orienta- 
tion cue is that as a surface is tilted away from the camera the edges in the direction of the 
tilt become foreshortened while the edges orthogonal to the tilt are not. This changes the 
edge orientation distribution and therefore the edge orientation distribution can be used as a 



cue for surface orientation. This effect is shown in figure 3.5 where the average HOG fea- 
ture is shown for various regions of tree trunk, forest floor, grass lawn, and patio tile. The 
cylindrical shape of the tree trunk is clearly indicated by the warping of the HOG feature 
as a function of position on the trunk. 

We consider the case of a surface with an isotropic edge distribution — the amount of 
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edge at any two edge orientations is equal. We do not expect this assumption to hold in 
any particular small surface patch, but we do expect a statistical relationship between HOG 
features and orientation. The assumption of isotropic surface textures provides a departure 
point for analysis. We are interested in the expected HOG feature for the edges projected 
onto the image plane when the surface is tilted relative to the image plane. 

We consider edges to be line segments on the image plane. Each edge has a length 
r and an orientation 0. The HOG feature is a histogram of edge orientations. We will 
normalize the HOG histogram so that it is a probability distribution on orientations. If P is 
a probability distribution (or probability density) on edges then the HOG feature H (0) is 
a probability distribution (or density) on 0. More specifically, H is a marginal of P onto 
in which edges are weighed by length. 

H(Q) = ^j P(0)E[r|0] (3.7) 

We consider a surface patch imaged by a perspective camera. A perspective camera induces 
the following map from three dimensional coordinates to image plane coordinates. 

x' = (fx)/z y' = (fy)/z (3.8) 

We assume a coordinate system on the surface patch such that each point on the surface 
patch has coordinates x s , y s . The image plane and surface coordinates can be selected so 
that we have the following map from surface coordinates to three dimensional coordinates 
where ^ is the angle between surface normal and the image plane normal. 

x = x s y = y s cos ^ z = zq + y s sin ^ (3.9) 
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Combining (|3.8[) and (|3.9[), and differentiating with respect to x s and y s gives the following 



with z = z$ + y s sin \l>. 

dx' = (tjdx s - y xs ^2 ) dy, 



(3.10) 
dy' (^)dy s -( fy«»*™* )dy 8 



z* 



We first consider the region near the center of the image. At the center of the image we 



have x = y = and (3.10) becomes the following. 

dx'=[-\dx s dy = ( J dy s (3.11) 

We now have the following map from the length R and orientation T on the surface to 
length r and orientation on the image plane. 

= tan - (sinrcos\l/, cosT) 



r 

z 



— V sin T cos 2 \l/ + cos 2 T 



For cos^I ^ Owe have a bijection between T and 0. Furthermore, for fixed and T we 
have a linear relationship between r and R. This gives the following. 



i/(0)cx Q(r(®))%viR\r(Q)}^ 



We assume that Q is isotropic, i.e., that Q(r)E [R\T] is a constant independent of T. We 
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then have the following. 

H(Q) oc — a/cos 2 * sin 2 T + cos 2 T 

(20 

dT 1 + tan 2 6 

dQ cos^fl + H^j 

V cos 2 * / 

We have no simple formula for if (O). However, we can see that if (6) is a continuous 
function of 6 satisfying the following. 



m = m = h{c£* 



H(n/2) = H(3n/2) = \ cos 2 * 

So we get the following where H m [ n is the minimum of H(Q) and if max is the maximum 
of if (9). 

ff m in/#max = COS 3 * (3.12) 



To justify the form of the orientation energy ( 3.6| ) we first note that in the same coordi 



nate system as (3.9) the equation for the surface plane can be written as follows. 

z = zq + y tan * (3.13) 

If we let b be the distance between the foci of the two cameras we have that the disparity d 
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equals bf/z. Multiplying (3.13) by bf /{zzq) gives the following. 



ZQ Z \ ZQ J \ Z 



d0 = ,,+ ,'^.y ( 3.i5) 



i = (fc _,'*^»y (3.16) 



B = -*^ (3.17) 

For the pixel p at the center of the image we have d(p) = do. We handle pixels outside 
of the center of the image by considering panning the camera to bring the desired point to 
the center and approximating panning the camera by translating the image. This gives the 
following general relation between the disparity plane parameter and the angle \& between 
the ray form the camera and the surface normal. 

B = _ d(p)t*a* (31g) 



In the orientation energy (3.6) we interpret /3g • H{p) as a predictor of — (tan\l/)// and 



we multiply by d(p) to get a predictor of B. We do not currently exploit (3.12). Doing so 
should result in a more refined texture cue. 

3.4 HOG Computation and Representation 

In order to use HOG features as a probability distribution on orientations, we want to 
compute HOG as a normalized histogram of edge orientations through the whole image. 
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Our HOG descriptor at each pixel is a 24-dimensional vector which describes the edge 
orientation distribution over the local image region around that pixel. (We compute a 8- 
dimensional HOG feature at each image scale, for 3 different scales to get a 24-dimensional 
feature vector at each pixel). The HOG computation includes the following 3 main steps: 

• Gradient Computation: 

The first step of calculation is the computation of the gradient values for all pixels 
in the left image. This is to simply apply the 1-D centered, point discrete deriva- 
tive mask in one of or both the horizontal and vertical directions: [—1,0, 1] and 
[-1,0, 1] T 

• Orientation Binning The second step of calculation involves creating the cell his- 
tograms. The cells themselves can either be square or radial in shape, and the his- 
togram channels are evenly spread over to 180 degrees or to 360 degrees, de- 
pending on whether the gradient is unsigned or signed. Here we use square cells, 
unsigned gradients in conjunction with 8 histogram channels. Each pixel within the 
cell casts a weighted vote for an orientation-based histogram channel based on the 
values found in the gradient computation. As for the vote weight, pixel contribution 
can either be the gradient magnitude itself, or some function of the magnitude; here 
we use the gradient magnitude. Specifically, we quantize the gradient orientation of 
pixels in a cell into 8 bins (corresponding to 8 equal ranges from to 180 degrees). 
In each bin we sum the magnitude of gradients in that bin. Then we normalize the 
magnitude of the resulting 8-dimensional vector of the bin. 

• Pyramidal HOG Descriptors We generate the HOG feature pyramid by computing 
HOG features of 3 different scales of the image pyramid (see Figure), corresponding 
to different sizes of the cell. The sizes of the cell at each scale are 8x8, 16x16 and 
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32x32 respectively. In this step we used the integral image trick for the purpose of 
computational efficiency. 
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Figure 3.5: HOG features for image regions. The amount of edge as a function of angle, a 
HOG feature, is averaged over different vertical and horizontal regions on various images. 
The different surface orientations of these regions affect the HOG features. We can see the 
cylindrical structure of tree trunks and the fact that the ground plane becomes more tilted 
as the distance increases. 
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Figure 3.6: We compute HOG features at each image scale for three different scales to 
obtain a 24-dimensional feature vector at each pixel. 



CHAPTER 4 
UNSUPERVISED CRF LEARNING FOR STEREO VISION WITH 

MONOCULAR CUES 



In this chapter, we apply the unsupervised CFR model we introduced in Section 1 .4 to the 
problem of unsupervised learning of a highly parameterized stereo vision model involving 
the monocular shape from texture cues - training the model parameters from unlabeled 



stereo pair training data( |Trinh and McAUester|2009 ). First, we describe in more details our 



PBP inference algorithm for this CRF stereo model. We then give a background review 
on contrastive divergence, a machine learning technique that plays a key component in our 
unsupervised CRF learning algorithm. Next, we present our unsupervised CRF learning of 
the slanted-plane stereo vision model. We conclude the chapter by reporting experimental 
results both for learning and inference on two different stereo datasets. 

4.1 Particle Belief Propagation for Plane Estimation 

Given a pair of images (X, Y), and a given segmentation of X, and a given setting of the 
model parameters, the inference problem is to find an assignment Z of plane parameters 



to superpixels so as to minimize the total energy in (3.2). The energy defines a Markov 
random field. More specifically, the match energy defines the unary potential on each 
superpixel independently and the smoothness energy defines the binary potential on pairs 
of adjacent superpixels. 



The complete inference algorithm is illustrated in Figure 4. 1 The steps of the algorithm 
include: 
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Initialization 



© 



© 



© 



Segment the left image 



Compute initial disp. Map 
for left and right images 



Occlusion detection by 
mutual consistency 



© 



Superpixel 

merging using 

plane parameters 



RANSAC 
plane fitting 
for initial z 



© 



Generate candidate 

values (particles) for 

each z. 



© 



Main loop 



Max-Product 
Particle BP to 
optimize E(z) 



8 



Compute HOG 
features for 
left image 



© 



© 



Update 



Return final z 



© 



Figure 4.1: Disparity plane inference algorithm using Particle-Belief Propagation. 



1. We compute a segmentation using the Felzenszwalb-Huttenlocher segmentation al- 



gorithm (Felzenszwalb and Huttenlocher 2004). Although this segmentation algo 



rithm has a parameter governing the number of superpixels, we do not attempt to 
tune this parameter with our training algorithm. A more principled approach would 
include the segmentation itself in the energy functional and tune the segmentation 
parameter along with the other parameters of the model. However by holding the 
segmentation parameter fixed, we can treat the segmentation as part of the input 
data. 



2. We run the Felzenszwalb-Huttenlocher's efficient loopy BP stereo algorithm ((Felzen 



szwalb and Huttenlocher 2006)) to compute the disparity maps for both images X 



andF. 



3. Mutual consistency check requires that for a particular pixel in the left image, the 
disparity values between the left and right disparity maps are consistent, i.e: 
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d L (p)=d R (p-d L (p)) (4.1) 



The pixel is considered occluded if (4.1 ) does not hold, and considered non-occluded 
otherwise. 

4. The plane fitting is performed in the disparity space, and is applied per superpixel 



to obtain an initial disparity plane. This is done robustly using RANSAC (Fischler 



and Bolles| 1981 ) on the disparity values of the non-occluded pixels only, (similar to 



dQ. Yang and Nistr|2008| )) 



5. This step is optional. In this step we implement a superpixel grouping algorithm 
similar to the one used in step 1. However each time we consider merging two 
adjacent superpixels, we also look at the difference in their disparity planes computed 
from step 4. Intuitively we prefer merging adjacent superpixels having very similar 
disparity planes. Note that if used, this step may change the initial segmentation 
generated by step 1 . Then the output is iteratively fed back into the plane fitting in 
step 4. 

6. At each pixel p in the left image X, we compute a HOG vector H(p) which is a 
24 dimensional feature vector consisting of three 8 dimensional normalized edge 
orientation histograms an 8 dimensional orientation histogram is computed at three 
different scales. 

7. Let Cj be a set of candidate values for z^ derived by repeatedly adding random noise 
to the current value of zp 

8. Run discrete max-product BP with the finite value set Q for each node i. 
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9. Set Z{ to be the best value for i found in step 8 and repeat. 

10. The iteration can be stopped after a fixed number of iterations or when the energy is 
no longer reduced. 

Note that the main loop of the algorithm, including steps 7, 8, 9, is very similar to 



the Particle -based Belief Propagation algorithm described in Section 2.4 which is closely 



related to previous work in (Ihler and McAllester 2009) and (Ko ller et al.||1999 ). This can 



be viewed as the max-product PBP algorithm. In Chapter [2[ we describe in details two 
other applications of the sum-product PBP inference algorithm to two inference problems 
in Computer Vision. 

4.2 Background Review on Contrastive Divergence 

4.2.1 Introduction to Contrastive Divergence 

Assume that our goal is to model the probability of the observed data x using a function of 
the form f(x, (3), where (3 is the model parameter vector. This probability can be defined 
as follows. 



P{x,0) = Jfl /(*,/*) ( 4 - 2 ) 



where Z(f3), commonly known as the partition function, is defined as 



Z(P) = / f(x,P)dx (4.3) 
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Given a training dataset X = (xi , X2, ...,Xtf), we want to learn the model parameter 
vector (3 so as to maximize the probability of the training data. 



N 1 
P(X,J3) = W^fi^P) (4-4) 



Taking the negative log of (4.4) will give us the following: 



1 N 
E(X,(3) = logoff) --£) log /(*<,£) (4.5) 



z=l 



We denote E(X, (3), which is the negative log of P(X, (3), to be the energy function. 



The partition function Z(f3) in equation (4.3) is, for many cases, a very complicated 



function which can be very hard to compute. In cases where this function is intractable, 
this leads to a situation where we are trying to minimize an energy function that we cannot 
evaluate. 

Contrastive divergence offers a solution to this problem. Briefly speaking, contrastive 
divergence is a tool that allows us to estimate the gradient of the energy function, even 
though we cannot evaluate the energy function itself. 

4.2.2 Formulation of Contrastive Divergence 

As previously mentioned, contrastive divergence provides us a way to estimate the gradient 
of our energy function E(X, (3) with respect to the model parameters (3, given the training 



dataset X. Taking the partial derivative of equation (4.5) with respect to (3 gives us the 
following derivation. 
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dE(X,/3) d log Z(P) 



dp 



dp 



1 N 

-Y 



i=l 



d log f(xj,P) 

dp 



(4.6) 



dhgZ(P) /dlogf(x,/3) 



dp 



dp 



(4.7) 



X 



where {.) x denotes the average over the training data X. In other words, it is the expec- 
tation of . given the training data distribution X. 



The first term on the right hand side of (4.6 ) involves the partition function Z{(5). Sub 



stituting equation (4.3 ) into this term, we have: 



dlog Z{P) 
dp 



1 dZ(/3) 
Z{p) dp 



(4.8) 



1 d 

~W)d~p 



f(x,P)dx 



(4.9) 



1 fdf(x,P) 



Z(P) J op 



dx 



(4.10) 



' W) 8 *'/^ 



Z(P) 



dp 



(4.11) 



/ 



P(x,P) Qj dx 



(4.12) 



dlog f{x,P) 
dp 



P(x,0) 



(4.13) 



From equation (4.13), it is not hard to see that, although the integration in (4.12) is 
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algebraically intractable, it can be numerically approximated by drawing samples from 
the proposal distribution P(x,(3). The next problem is that the distribution P(x,(3) is 
actually unknown, as we do not know the value of the partition function. In order to resolve 
this, we can use a Markov Chain Monte Carlo (MCMC) sampling method to simulate 
draws from the proposal distribution by drawing directly from the target distribution (the 
training data distribution), then transforming these samples to the proposal distribution. 
The transformation is done by generating a Markov chain starting from the drawn sample 
at the target distribution to the proposal distribution. Please recall that in the Markov chain, 
each state x t+l only depends on the previous state x t . At each time step t, a new value x' 
is accepted for the next state x t+l if a random number drawn from £7(0, 1) satisfies: 



. \ P{x' ,P)Q{x\x') 
a ^{Pix^^Qix^xt)' 

P(x' 0) 
This involves computing the ratio , v t ' , which is possible since the partition function 



cancels out. The approximated version of equation (4.6) can be written as follows 



dE(X,P) /dlogf(x,P)\ /dlogf(x,/3)\ (M1) 



df3 \ d/3 I X n \ dp /x 

where X n represents the set of samples drawn from the training data after n steps of 
MCMC. 

The use of MCMC sampling strategy, however, can be very inefficient and not prac- 
tical, since it may take a very long time for the Markov chain to mix and approximately 



converges. Hinton in (Hinton 2002) asserts that only a few MCMC steps would be sufficient 
to calculate an approximate gradient. The intuition behind this is that after a few iterations 
the data will have moved from the target distribution (i.e. that of the training data) towards 
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the proposed distribution, and so give an idea in which direction the proposed distribution 
should move to better model the training data. Empirically, Hinton has found that even 1 
step of MCMC is often sufficient for the algorithm to converge. 

That said, as we use gradient descent in order to minimize our energy function, the 
contrastive divergence parameter update equation should look like the following: 



where 7 is the step size value which is allowed to change at every iteration, based on 
convergence time and stability. 



4.3 Our Unsupervised CRF Learning Algorithm 



As described earlier in Section 3.1 our stereo model involves three terms — a correspon- 
dence energy measuring the degree to which the left and right images agree under the in- 
duced disparity map, a smoothness energy measuring the smoothness of the induced depth 
map, and an texture energy measuring the degree to which the surface orientation at each 
point agrees with a certain (monocular) texture based surface orientation cue. For surface 



orientation cue we use histograms of oriented gradients (HOG) ( Dalai and Triggs|[2005a ) 



The stereo model itself and the monocular surface cues (as part of the model) are both 
viewed as the objects being trained. Our energy functional includes in total 62 parame- 
ters - 10 correspondence parameters, 2 smoothness parameters, and 50 texture parameters. 
Among these parameters, we train 50 of them are trained using gradient descent. The 
gradient of the energy function w.r.t to each parameter is obtained using contrastive diver- 
gence, the method described in the previous section. The other 12 parameters, including 
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10 matching parameters (A^ in (3.5)) and 2 texture parameters (fij^ and /3g in (3.6)), are 
computed by least squares regression. 

Our approach to unsupervised learning is based on maximizing conditional likelihood. 
In particular, we consider a general conditional probability model Pg(u\x) over arbitrary 
variables x and u and defined in terms of a parameter vector (3. Analogously, CRFs model 



( |Lafferty et al.||200~T ) defines a conditional probability of the dependent variables u given 



the exogenous variables x and (importantly) does not model the distribution of the exoge- 
nous variables. In the case of stereo vision one might expect that it is easier to model the 
probability distribution of the right image given the left image than to model a probability 
distribution over images. 



In a closely related earlier work by Zhang and Seitz (Zhang and Seitz 2007 ), the authors 
also used a similar CRF model for the classical pixel-based stereo algorithm. Although the 
training was also based on maximizing conditional likelihood, they worked with a much 
lower-dimensional model (5 parameters) and tuned these parameters separately to each 



single stereo pair as in (4.16). 



Pf = argmax Pg{ui\xi) (4.16) 

The five parameters are tuned to each individual input stereo pair, although the method 
could be used to tune a single parameter setting over a corpus of stereo pairs. The main 
difference between their work and ours is that we train highly parameterized monocular 
depth cues. Another difference is that we formulate a general CRF- like model for unsu- 
pervised learning based on maximizing conditional likelihood and avoid the need for the 
independence assumptions used by Zhang and Seitz by using contrastive divergence — a 
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general method for optimizing loopy CRFs (Hinton 20021 |Carreira-Perpinan and Hinton 



2005 1. 



There is also related work on learning highly parameterized monocular depth cues 
by Saxena et al. in (Ashutosh Saxena and Ng 2007b|a ), as well as Hoiem et al. in 
( |Derek Hoiem|[2005] ). The main difference between these methods and ours is that we 
use unsupervised learning while they use ground truth data to train their system. In Chap- 
ter [• 



4.4 we will show that our algorithm with unsupervised training outperformed the results 



in (Ashutosh Saxena and Ng 2007 a| ). One might argue that stereo pairs constitute super- 
vised training of monocular depth cues. A standard stereo depth algorithm could be used to 
infer a depth map for each pair which could then be used in a supervised learning mode to 
train monocular depth cues. However, we demonstrate that training monocular depth cues 
from stereo pair data improves stereo depth estimation. Hence the method can be legiti- 
mately viewed as unsupervised learning of a stereo depth. Also the general formulation of 
unsupervised learning by maximizing conditional likelihood, like the shift from MRFs to 
CRFs, may have significance beyond computer vision. 

Other related work includes that of Scharstein and Pal ( |Scharstein and PaT||2007| ) and 
Kong and Tao ( Kong and Tao|[2004 ). In these cases somewhat more highly parameterized 
stereo models are trained using methods developed for general CRFs. However, the training 
uses ground truth depth data rather than unlabeled stereo pairs. 

We learn MRF parameters using contrastive divergence ( |Hinton 2002 ; Carreira-Perpinan 
and Hinton|2005 ), a general MRF learning algorithm capable of training large models. An- 



other work that also described learning on MRF using contrastive divergence is the Field 
of Experts system by Roth and Black in ( |Roth and Black|2005| ). In ( |Roth and Black|2005| ), 
the authors describe an extension of the traditional MRF model by learning potential func- 
tions over extended pixel neighborhoods, which they then implemented with two vision 
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applications: image denoising and image impainting. To training the model parameters, 
they also used contrastive divergence update, aiming at minimizing the KL-divergence be- 
tween the model distribution and the data distribution. The MCMC sampling strategy was 
used to approximate the expectation over model distribution. The data distribution is easy 
to compute by computing the average over training data. Similar to other aforementioned 
related work, they also used training data with ground truth label. 



Our training approach is actually more related to the hidden CRF model ( Quattoni et al. 



2007 1 ) than the standard CRF. Our model includes one more variable, the latent variable y 
which is not observed in the training data. Specifically in our slanted plane stereo model, 
x denoted segmented left image, u denoted a right image, and y was an assignment of a 
plane to each superpixel of x. y is not in the training data, but it is part of the model. 
However, in this chapter we aim at describing our model in a more general level rather 



than a model specifically applied for stereo vision, as defined by (4.17). The conditional 
probability model Pg{u\x) now is defined not only over x and u but also by an arbitrary 
latent variable y. 

Pp{u\x)=Y,Pfi(u>v\x) ( 4 - 17 ) 

y 

Given a set of stereo image pairs as training data (xi , u\), . . . (xjy, u^-) conditional 
EM is an algorithm for locally optimizing the parameter vector (3 so as to maximize the 
probability of the u values given the x values in the training data, as formulated in equation 
( |4 . 1 8 > (This is exactly equation (1.16) from Section 1.4). 



TV 
(3* = argmax V^ ha. Pa (ui\xi) (4.18) 

Conditional EM is a straightforward modification of EM and is defined by the following 



two updates where (3 is initialized with domain specific heuristics. 
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p i(y)-= p p(y\xi; 



U4 



(4.19) 



TV 



p:= argmax J^E y ^p. [In Pp(iii,yi, 

P i=l 






(4.20) 



Update ( |4.19| ) is called the E step and update ( |4.20| ) is called the M step. Hard EM, 
also known as Viterbi training, works with the single most likely (hard) value of y rather 



than the (soft) distribution P{ defined by (4.19 ). Hard conditional EM locally optimizes the 



following version of (4. 1 8 ). 



N 



P* 



argmax > max In Pq {u\ , y \xi) 



(4.21) 



i=l 



Hard conditional EM is defined to be the process of iterating the updates (4.22) and 
(4.23 ) below which can be interpreted as hard versions of (4.19 ) and ( 4.20[ ). 



Ui := argmax Pp(ui, y 



X } 



(4.22) 



N 



/3 := argmax J^lnPg^,^ 

P i=l 



X j 



(4.23) 



We will call (422) the hard E step and (423) the hard M step. Updates (422) and 



(4.23 ) are both coordinate ascent steps for the objective defined by (4.21 ). However, we 



refer to (4.22) and (4.23) as hard conditional EM rather than simply "coordinate ascent" 



because of the clear analogy between (4.21 ), (4.22), (4.23) and (4.18), (4.19), (4.20) 
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The parameter vector (5 is composed of two components (3 = ([3 y , (3 U ) where /3y 
parameterizes P(y\x) and f3 u parameterizes P(u\x, y). Note that the probability model on 



the right hand side of (4.22 ) can be further expanded as follows: 



V, 



argmax Pp y (y\x t ) Pp u (ui\x i: y) 



argmaxln [Pp^x^Pp^u^x^y) 

argmin Er (x t , y) + Ep u {x{ ,u i: y) + In Zp u (x, y) + In Za (x) (4.24) 



We have that: 



Z (3y( X ) 



y 



(-Ey(x,y,l3y) ) 



(4.25) 



z p u ( x >y) = J2 e 



(-E u (x,u,y,/3u) ) 



(4.26) 



where equation (4.25) corresponds to the sum of (3.4) and (3.6), equation (4.26) corre 



sponds to (3.5). Since (4.25) does not depend on y, we can eliminate it from (4.24). It is 



less obvious to show that (4.26) does not depend on y. 



z f3 u ( x >y) = Yl e 



{-E u (x,u,y,f3 u ) ) 



J2 e -&pex($ x (p)-$ u (p+d(p,y)) ) 2 ) 



(4.27) 
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We assume that the mapping from x to u in (|4.27|) is a bijection, i.e. each pixel p in 



x only maps to one unique pixel in u (since x and u have the same size, this condition is 



sufficient to guarantee bijection). We can consider the matching cost from equation (3.5 ) 
as a distance function between the left image x and one of its prediction x, where x is a 



function of a pair (u, y) . Note that (4.26 ) is a sum over all possible values of u, and pairing 



all possible u with y can generate all possible values of x. We can rewrite (4.26) as follows: 



z hS?>V) = Yl e( 

u 

= E 



-E u (x,u,y,fiu) ) 



i(- E x( x $>Pu) 



(4.28) 



It follows that (4.26) does not depend on y either, and therefore can also be eliminated 



from (4.24). We can rewrite (4.24) exactly as: 



y { := aTgmmE /3 (x i ,y) + E (3u (x l ,u l ,y) 
V 



(4.29) 



Each of the terms in the right hand side of (4.29) is a CRF. In the case of the slanted 



plane stereo model, the hard E step (4.24 ) is implemented using a stereo inference algorithm 



which computes yi by minimizing the corresponding energy functionals. The inference 



algorithm is described in Section 4.1 



Our implementation of the hard M step relies on a factorization of the probability model 
into two conditional probability models each of which is defined by an energy functional. 
Unlike CRFs, we do not require the energy functional to be linear in the model parameters. 
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D pu,p y ( u ^\ x ) = P l3 y (y\ x ) p i3 u (u\x,y) (4.30) 



J-E y (x,y,(3y) ) 



Z y (a;,/3y) = ^ e (-Ey(x,y,p y ) ) 
y 

(-E u (x,u,y,p u ) ) 
^/3>M) = —t7 t^— (4.32) 



Z u (x,y,p u ) = Yl e 



{-E u (x,u,y,P u ) ) 



Given this factorization of the model, the hard M step (4.23) can be written as the 
following pair of updates. 



P y := argmaxJJP^(j/i|xi) 

h i 

: = axgmax^lnP^^I^) 

h i 

:= argminy^ Ey(xi,yi,p y ) 

h i 

/3 U := argmaxY[ Pi3 u (ui\xi,yi) 
Pu i 

:= argm&x^ln Pp u (ui\xi,yi) 

Pu i 

:= axgminy^ E u (xi,Ui,yi,/3 u ) (4.33) 

Pu i 
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Let L abbreviate the quantity being maximized in the right hand side of (0331) and let 
Ei(y) abbreviate E y {xi, yi, (5). We can express the gradient of L as follows. 



N 



W Py L = Yl { E y~Py(y\xi,P) 
i=l 



Vp^y) -Vp y Ei( yi 



(4.34) 



A similar equation holds for (4.33 ). 



The first term in the right hand side of (6.22 ) is the expectation over the model distribu- 
tion, while the second term is the expectation over the training data. The first term involves 
Pa (y\xi), which is a very complicated probability distribution and usually is extremely 
difficult to compute. However in such cases, we can estimate the expectation of such distri- 
butions by taking the average over a sufficient amount of samples, obtained by sampling y 
from Pp (y\xi) using an MCMC sampling process. The regular MCMC sampling method 
is performed by running a sequence of Metropolis or Metropolis-Hastings sampling steps 
until the Markov chain mixes - i.e. reaching the stationary distribution. Usually it is a 
difficult problem to determine how many steps are needed to converge to the stationary dis- 
tribution within an acceptable error. Contrastive divergence proposed by G. Hinton et. al. 
(Hinton 2002; Carreira-Perpinan and Hinton 2005) is a technique that allows us to sample 
z using only one or two MCMC steps rather than a long running MCMC process. Since we 



can estimate V« L, we can then optimize (4.33 ), and similarly (4.33 ), by gradient descent 



For the experiments reported here we use contrastive divergence ( Hinton|2002 Carreira- 
Perpinan and Hinton|2005[ ) to sample y rather than a long running MCMC process. In con- 
trastive divergence we initialize y to be y^ and then perform only a few (one or two) MCMC 
updates to get a sample of y. Contrastive divergence can be motivated by the observation 
that if yi is assumed to be drawn at random from Pg(y\xi) then the expected contrastive 
divergence update is zero. So as j3 better fits the pairs (xi, yi) one expects the contrastive 
divergence gradient estimate to tend to zero. Furthermore, because only a few updates are 
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used in the MCMC process, contrastive divergence runs faster and with lower variance than 
a longer running MCMC process. Contrastive divergence yields a consistent estimator for 
(3 whenever the expected update direction is the gradient of a convex potential function. 

We will demonstrate the experimental results in the next chapter to show that training 
monocular cues from stereo pair data alone improves stereo depth estimation. Our unsuper- 
vised learning algorithm implicitly trains shape from texture monocular surface orientation 
cues. Moreover, our stereo model with texture cues, only by unsupervised training, out- 
performed the results in ( Ashutosh Saxena and Ng|2007a ). Throughout this chapter, stereo 



vision was considered as a simple setting to investigate unsupervised learning and hence 
seems a good place to start. However we believe that our approach to unsupervised learn- 
ing based on maximizing conditional likelihood can be generalized to other, maybe much 
more sophisticated models. 

4.4 Experimental Results 

We implement two training methods in our experiments — supervised and unsupervised. 
For each of the supervised and unsupervised training methods we train both a version with 
texture cues for surface orientation and a version without such cues. In supervised training 
we set z\ (for each training image) by fitting a plane in each segment to the ground truth 
disparities for that segment. We then train the model using a contrastive divergence imple- 



mentation of the hard M step (4.23 ) which we describe in more detail below. In supervised 



training we use only a single setting of z and run one iteration of (4.23). In unsupervised 



training we use the same separation into training and test pairs but do not use ground truth 



on the training pairs. Instead we iterate (4.22) and (4.23) six times starting with initial 
values for the parameters. 



77 



We use the inference algorithm described in section |4~T| to implement the hard E step 



(4.22). This uses a form of max-product particle belief propagation. Given an assignment 
z of a plane to each segment, we propose 15 additional candidate planes for each segment 
by adding Gaussian noise to the plane specified by z. The plane parameters A and B have 
units of pixels of disparity per pixel in the image, and hence are dimensionless. Typical 
values of | A | and \B\ are from .1 to 1. In the proposal distribution we use Gaussian noise 
with a standard deviation of .007 for each of A and B and use a deviation of . 1 pixels for 
C. We perform six rounds of proposing and selecting. 

We implement the hard M step ( 4.23[ ) by first breaking it down into (4.33) for training 
P(z\x, (3 Z ) and (4.33 ) for training P(y\x, z, (3 y ). The form of the match energy (a sim- 



ple quadratic energy) allows a closed form solution for (4.33). We implement ( 4.33[ ) by 
gradient descent using a contrastive divergence approximation of the gradient in ( 6.22[ ). 
We perform 8 gradient descent parameter updates with a constant learning rate. To esti- 



mate the expectation in (6.22) in each parameter update we generate 10 alternative plane 
assignments using single MCMC stochastic step starting at z and accepting or rejecting 
Gaussian noise added once to each plane. The MCMC process proposes a new plane for 
each segment by adding Gaussian noise and then accepting or rejecting that proposal using 
the standard Metropolis rejection rule. 



4.5 Experimental Results with the Middlebury Dataset 
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all 
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2.58 


4.66 


11.8 


0.47 


0.64 


6.1 


6.72 


6.98 


16.1 


6.93 


9.33 


16.6 


7.41% 



Table 4.1: Performance on the Middlebury stereo evaluation. The numbers shown are for 
unsupervised training with texture features. 
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Table |4.1| shows the performance of our system on the Middlebury stereo evaluation 



(version 2). The numbers shown are for unsupervised training with texture features. In 
this case all four images where used as unsupervised training data (ground truth disparities 



were not used in training). Figure 4.2 shows the inferred depth maps for the Middlebury 
images at various points in the parameter training. The figure shows a clear improvement 
as the parameters are trained. 




Left image Iteration 1 Iteration 3 Iteration 5 

Figure 4.2: Improvement with training on the Middlebury dataset. 
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RMS 

Disparity Error 
(pixels) 


Average 
Error 

1 !ogio z ~ logio Z\ 


Saxena et al. (Ashutosh Saxena and Ng 2007a) 




.074 


Unsuper., Notexture 


1.286 


.075 


Unsuper., Texture 


1.1608 


.0723 


Super., Notexture 


1.1142 


.0709 


Super., Texture 


1.0319 


.0686 



Table 4.2: RMS disparity error (in pixels) and average error (average base 10 logarithm 
of the multiplicative error) on the Stanford stereo pairs for four versions of our systems 
plus the best reported result from ( |Ashutosh Saxena and Ng|[2007a| ) on this data. Each 
system was either trained using the ground truth depth map (supervised) or trained purely 
from unlabeled stereo pairs (unsupervised) and either used texture cues (Texture) for sur- 
face orientation or did not (Notexture). Note that texture information helps improve the 
performance in both supervised and unsupervised cases. 



4.6 Experimental Results with the Stanford Dataset 



We have also run experiments on a set of stereo image pairs taken from the Stanford color 
stereo dataset [J] which has been used to train monocular depth estimation (Ashutosh Saxena 
2003 1 Ashutosh Saxena and Ng 2007b|a ). The images cover different types of outdoor 
scenes (buildings, grass, forests, trees, bushes, etc.) and some indoor scenes. 

First we performed epipolar rectification on this dataset. Given a pair of stereo im- 
ages, rectification computes a transformation matrix for each image plane such that pairs 
of conjugate epipolar lines become collinear and parallel to one of the image axes (usually 
the horizontal one). The rectified images can be thought of as acquired by a new stereo 
rig, obtained by rotating the original cameras. The important advantage of rectification is 
that computing stereo correspondences is made simpler, because search is done along the 
horizontal lines of the rectified images. 

Standard rectification methods such as ( Irsara and Fusiello|2006[ ) rely on detecting and 



1. http://ai.stanford.edu/ asaxena/learningdepth/data 
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matching a sparse set of feature points to compute the epipolar geometry between the left 
and right image. The rectification then is achieved usually by minimizing a measure of 
distortion - a score function involving only the point correspondences. 

Using a rectification kit from Fusiello et al. ( Irsara and Fusiello||2006 ) failed in 57 out 



of 257 stereo pairs in the Stanford dataset (The failed cases are cases where the rectifica- 
tion score function exceeded a specified threshold). Therefore we came up with a different 
rectification approach that we call '"dense rectification'". First we observe that all stereo 
pairs in the dataset were taken from a stereo rig with almost the same calibration parame- 
ters. We then aim at determining one or only a few rectification solutions for all pairs in the 
dataset, rather than computing a separate rectification for each pair. Besides, since epipolar 
rectification is considered an important preprocessing stage of dense stereo matching, the 
ultimate goal of rectification is also to optimize the dense stereo matching score. Therefore 
our method strives to find the rectification solutions that directly serve for this purpose: 
minimizing the dense stereo energy score rather than the feature point correspondence dis- 
tortion score. To compute the dense stereo energy score for each pair, we simply use the 
efficient stereo algorithm in ( Felzenszwalb and Huttenlocher||2006 ). Our approach can be 



summarized as follows: 



Run rectification code from Fusiello et al. (Irs ara and FusieUo||2006| ) on the whole 
dataset. 

Use stereo energy to pick the best 200 image pairs. (Consider the other 57 images 
outliers) 

Do K-means clustering on the 200 transformation matrices of these 200 image pairs 
with (K = 2,3,4). 

Pick best K = 2. Compute the 2 mean transformation matrices. 
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• Apply each of the transformation matrices to rectify the 257 images. For each image 
we pick the rectified one with better stereo energy score (instead of the rectification 
score). 

• We obtain 257 rectified image pairs. 

With our rectification algorithm, we successfully rectified all stereo image pairs in the 
dataset. The rectified dataset has been put online for public use^] Our results are directly 
comparable to ( |Ashutosh Saxena and Ng||2007a[ ), i.e we use the same training set and test 



set (193 pairs for training, 64 pairs for testing). Table |4.2| shows that we outperform their 
results. 

4.7 Conclusion 

In many applications we would like to be able to build systems that learn from data col- 
lected from mechanical devices such as microphones and cameras. Stereo vision provides 
perhaps the simplest setting in which to study unsupervised learning. We have formu- 
lated an approach to unsupervised learning based on maximizing conditional likelihood 
and demonstrated its use for unsupervised learning of stereo depth with monocular depth 
cues. Ultimately we are interested in learning highly parameterized sophisticated mod- 
els including, perhaps, models of surface types, shape from shading, albedo smoothness 
priors, lighting smoothness priors, and even object pose models. We believe that unsuper- 
vised learning based on maximizing conditional likelihood can be scaled to much more 
sophisticated models than those demonstrated here. 



2. http://ttic.uchicago.edu/ ntrinh/shared/stanford/ 
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Figure 4.3: Several depth map results on the rectified Stanford stereo dataset. The first 2 
row are some example images. The second 2 rows are our output depth maps. 



CHAPTER 5 
UNSUPERVISED CRF LEARNING FOR MONOCULAR DEPTH 

ESTIMATION 

In this chapter we present our unsupervised CRF learning approach to the problem of 
monocular depth estimation (MDE). MDE is the problem of predicting the depth at each 
pixel of a single input image. This is a much more challenging problem than stereo vision. 
There exist an intrinsic ambiguity between local image information and the 3-D location 
of the real- world object, due to perspective projection. In other words, a point in the image 
plane can correspond to any 3-D point lying on a line connecting the camera center and 
that image point. 



In recent work by Saxena et al. (Ashutosh Saxena 2005; Ashutosh Saxena and Ng 



2007b |a ) an MDE system is trained from images paired with laser range finder depth maps. 
Here we are interested in training an MDE system using only stereo pairs — pairs of images 
taken from slightly different angles. Stereo pair training for MDE is interesting for a variety 
of reasons. First, in some cases it may be easier to obtain stereo pairs than to obtain images 
paired with laser range finder data. Second, learning MDE from stereo pairs can be viewed 
as a first step toward learning MDE from video sequences. Video is easily obtained and 
widely available and could provide an enormous volume of training data. Third, stereo pair 
training for MDE serves as a case study in unsupervised two-view learning which might 
provide a model for this form of learning in other applications. 

We train a monocular depth estimator using the same stereo pair database as used in 



(Ashutosh Saxena and Ng 2007a) but without the use of ground truth depth maps. For 
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learning we use basically the same unsupervised CRF learning algorithm we applied for 
Stereo Vision with Monocular Cues in Chapter |4} Our hard EM algorithm first estimates 
depth using stereo alone (E step), then trains the monocular predictor (M step), then re- 
estimates depth using both stereo and monocular cues, then retrains the monocular pre- 
dictor, and so on. While most of the performance of the monocular predictor is achieved 
by training on the initial stereo depth estimates, a modest improvement is seen for both 
view prediction and for ground truth prediction at later EM iterations. The performance on 
ground-truth prediction is similar to that reported in ( Ashutosh Saxena and Ng||2007a ) in 
spite of the absence of ground truth in training. 

In this chapter, first we present our depth estimation model with monocular depth cues. 
Next, we introduce the notion of view prediction and its relation to the conditional likeli- 
hood that the learning algorithm wants to maximize. We then describe the application of 
our unsupervised CRF learning algorithm for this specific model, aiming at training the 
MDE. The last section demonstrates experimental results. 

5.1 The model 

Our model is simply an extension of the standard dense stereo model. 



E w ^\l^,d) 



^dE P ^(\\Y l il \p)-Yp\p + d l (p))\\ 2 ,r d ) 
+ x sT,p, q eN( P ) min(\d(p) -d(q)\,T s ) (5-1) 

k +J2 P (d(p)-w-X( l )(p)) 2 



where (J*- 1 -*, 1^ >) is a stereo pair, d is the depth map of 1^ >, d(p) is the depth value 
assigned to pixel p, Y^ (p) is the image feature vector corresponding to pixel p in I'-v, 
jW(p) is the monocular feature vector of p, N(p) denotes the set of neighbors of pixel 
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p, and w is the MDE parameter. For this current model, we assume that other parameters 
Xfl, A s , Td, r^ are fixed, and the only parameter we want to learn is w. 



Note that (5.1) defines an MRF, in which the features are stereo, smoothness and 



monocular features, the hidden labels are the pixel depth values. The first term in (5.1 ) 
is the match energy measuring how well the disparity assignment d agrees with the input 
image pair, the second term is the smoothness energy enforcing the smoothness assump- 
tion, and the third term is the MDE term. For MDE we want to construct a disparity map d 
given only the first view l' K We assume a feature vector X' '(p) at each pixel p of /' ). 



A particular choice of features will be discussed in Section 5.4 but we can think of X' > (p) 
as containing a constant feature (a bias feature), information about the y component of the 
pixel p, and the value of various filters applied to image l' 1 -* at the pixel p. Our MDE 
energy term now is the sum squared difference between the assigned disparity d(p) and a 
locally predicted disparity w ■ X(w). 

We now define a monocular disparity predictor rfj^(/' ') determined by a parameters 
vector w as follows. 



4//W) = argmin £(d(p)-wlW(p)) 2 + A s ^ min(\d(p) - d{q)\ ,r s 
d p p,q&N( P ) 



(5.2) 



The second term here is exactly the robust L\ smoothing term in (5.1 ). (5.2) itself also 
defines an MRF, without the stereo feature. 
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5.2 Stereo Pair View Prediction and Conditional Likelihood 

Here we introduce the notion of view prediction. In view prediction the goal is to predict the 
second view given the first. View prediction training can be expressed with the following 
schema where the distortion and complexity functions are application specific. 

M 
/* = argmin \. Distortion I f lx^ J , x\ J + AComplexity(/) (5.3) 



In the system developed here we take / f x^> ) to be an image derived by applying 
an inferred disparity map to the image x^ 1 ' and the distortion function is simply the pixel- 
wise squared error (the L^ distance) between a;' > and / f x^ > J . A log-loss version of view 
prediction can be formulated as follows where P w {x^ '\x^ >) is a conditional probability 
(or probability density) parameterized by w. 

M 1 
w* = argmin y. m 7 c + AComplexity(w) (5.4) 

w f~i Pw (xW\x( l n 

For stereo pairs we might define P w f a?' ' |x' > J by inferring a disparity map from x^ ' 
under parameter setting w and then using the inferred disparity map to predict a probability 
density in color space at each pixel of x^ > . Note that view prediction training as defined 
by (5.3) or (5.4) is meaningful without labels. 

We consider stereo pairs of the form (/' > , l' >) where P > and j' 2 ' are images each 
of which is an assignment of a color vector to pixels. We write l(i)> c (p) for the color c 

nponent at pixel p in image I^> . Let d denote a disparity map — a function that assigns 



component at pixel p in image I^> . Let d denote a disparity map — a function that assi,_ 
a disparity d(p) to each pixel p. Disparity defines a mapping between the two images. The 
position (x, y) in P l > corresponds to the position (x, y + d(x, y)) in image /' <*. We let p 
range over pixels where each pixel is defined by a pair of integer coordinates (x, y) and we 
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write p + d(p) as an abbreviation for (x, y + d(x, y)) where (x, y) is the pixel p. We assume 
that d(p) is integral. We will mostly work with disparity d rather than depth Z where d and 
Z are related by Z = bf/d where b is the baseline (the distance between the two cameras) 
and / is the focal length. We assume that b and / are known. 

Given a disparity map d and a reference image I^> we construct a prediction P > which 
approximately satisfies the following equation. 



(2) c 

I^ c (p + d{p)) = ^(/ (1) ' c (p) - / (1) ' c ) + J (2) ' c (5.5) 

Here JW' C and en > iC are the mean and standard deviation of all the color c pixel values 
in image i' 1 -*. The quantities Jv 2 )> c and a^ >' c are defined similarly for image /' > . The 
use of pixel means and variances for the two images corrects for different biases and gains, 
in each color channel separately, of the two cameras. We could have assumed knowledge 
of camera bias and gain of the two cameras rather than knowledge of pixel mean and 
variance of l' > . But in any case the mean and variance of the pixels in l( > is only a very 



small amount of information about /' ' . Note that equation (5.5) both overspecifies and 
underspecifies Jv 2 J as a function of 1^ > and d. Overdetermined values in Jv 2 ) are set to the 
pixel value in /W with largest disparity (the closest point) and missing values in 1^ > are 
filled in from neighboring values preferring values derived from pixels with low disparity 
(farthest points). 

The overall view predictor f w (/' M is then defined by (5.5) applied to I^ 1 ' and 
dw M )• We will work with the following distortion function where P is the number 
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of pixels and C is the number of colors. 

Distortion^, ,<*>) = JL E d^(p)-I^(p)f (56) 

p,c (a^)> c ) 

Here we have defined distortion to be a percentage of variance. Note that if we take 
1^ > to be simply the mean pixel value of /' ' then we get exactly 100% distortion. 

We have pointed out that the view prediction error ( 5.6[ ) corresponds directly to the 



conditional probability P w (l^)\I^> ], i.e. the probability of the right image given the 
left image. This probability is exactly the conditional probability in our unsupervised CRF 
model: 

N 
(5* = argmax \. ^ n P/3( u i\ x i) (5-7) 

I 3 i=\ 

where u corresponds to the right image /' \ x corresponds to the left image I^ 1 ), and 
(3 is the model parameter, which corresponds to w. 

5.3 Unsupervised CRF Learning for MDE Model 

We are now interested in training the parameter vector w so as to minimize the expected 



distortion (5.6 ) over fresh stereo pairs. We assume a collection of training pairs 



(1) (2)\ (M) J2) 



We now formulate the hard EM algorithm which first estimates disparity using a clas- 
sical (untrained) stereo algorithm (E step), then trains a monocular estimator from the in- 
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ferred disparity map (M step), then re-estimates disparity using both binocular and monoc- 
ular cues, retrains the monocular predictor, and so on. We formulate this EM bootstrap 
algorithm as coordinate descent on a joint energy function of both the disparity map d and 
the parameters w. Formulating this algorithm as coordinate descent guarantees conver- 
gence. 



The joint monocular and binocular energy function was defined in (5.1 ). In our exper 



iments we take Y- (p) to consist of the color vector at p and the gradient of each color 
channel at p for a total of 3 + 6 = 9 feature values for three colors |^| The hard EM learning 
algorithm is then defined by the following optimization problem. 



w 



N 

argmin \. min£J u; (/- , J. , d) (5.8) 

w ^~ d 



Our training algorithm is coordinate descent on (5.8). We alternatively optimize the 



"coordinates" w and rfj. We initialize d{ to be the disparity map minimizing the first two 



terms of (5. 1 ). This corresponds to classical (untrained) binocular inference of the disparity 



map. Given an initial value for the disparity maps d{ we alternate the following two updates. 



TV 
w : = argmin J2 E w (lW , I& , d t ) (5.9) 

w i=l 
N 
= argmin V V (^ (p) - w ■ X\ (p) ) 2 
w i=l p 

d % := argmin E w (if \lf\d) (5.10) 



1. Before computing color vector and gradient features the images are normalized to remove 
bias and gain effects by subtracting the mean color and then dividing each color channel by the 
standard deviation of that color channel. 
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Table 5.1: Quantitative comparison between the ground plane baseline and our model using 
the average testing error per pixel as functions of the number of iterations. For each model 
we report: (a) The average view prediction error, measuring the predicted right frame and 
the ground truth right frame, i.e. the quantity on the left side of ( 5.6[ ). (b) The RMS disparity 
error compared with the ground truth (in pixels), measuring the average difference per pixel 
in disparity value between the predicted disparity and the ground truth disparity. 



Note that (5.9) is a least squares regression problem that can be solved exactly with 



standard codes. Also note that regularization of w can easily be incorporated into (|5.9[). 



We approximately solve (5.10) using the efficient loopy BP algorithm of Felzenszwalb 



and Huttenlocher (Felzenszwalb and Huttenlocher 2006). To relate this to the learning 
algorithm described earlier in Chapter|4} please note that equation ( 5.9 ) corresponds exactly 
to equation (4.23), and (5.10) corresponds exactly to equation ( 4.22| ). 

For the experiments reported in SectiorJ54| parameters of the algorithm other than w 
are set by hand. 



5.4 Experimental Results 

We ran the EM bootstrap training algorithm on a set of rectified stereo pairs taken from a 
color stereo datasetjjThe images cover different types of outdoor scenes (buildings, grass, 
forests, trees, bushes, etc.), and some indoor scenes. The images were epipolar rectified 



using a rectification kit from Fusiello et al. (Irsara and Fusiello 2006 ). Initial disparity maps 
were computed for all images using Felzenszwalb and Huttenlocher' s efficient loopy BP 



2. The dataset is available at http://ai.stanford.edu/ asaxena/learningdepth/data. 
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Left image Stereo depthmaps Monocular depthmaps Mono+Stereo depthmaps 

Figure 5.1: Resulting disparity maps from (a) the initial stereo algorithm, (b) the trained 
monocular predictor, (c) the combined monocular and binocular predictor. 



algorithm on the first two energy terms of ( |5.1[ ). By examining the resulting rectifications 
and depth maps it was clear that some rectifications had failed. At this point we removed 
from the dataset all pairs for which the energy value achieved by loopy BP was above a 
specified threshold. This left 204 out of an original 250 stereo pairs on which to train the 
monocular depth estimator. 

The monocular feature vector X' '{p) consists of local image features computed at 3 



different scales, which capture the texture and color cues (similar to ( Ashutosh Saxena and 



Ng~||2007a )). We also used the pixel image coordinates as features and included 1 constant 
feature for a bias term. We have X^>{p) e B? (17 image features at each scale for 3 



scales, plus 2 spatial features and 1 bias feature). Following (Ashutosh Saxena and Ng 



2007 a| ) we made disjoint copies of this feature set for different height levels in the image so 
that the monocular predictor is trained separately for each height level. We use 40 discrete 
height levels. 
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Training the monocular predictor separately at different height levels exploits the fact 
that these images were taken from a single robot with a fairly stable ground plain orientation 
relative to the cameras. The stability of the ground plane in these images suggests a baseline 
monocular prediction in which the disparity at each height level is simply predicted to be 
the average disparity for that height level across all training data as estimated in the learning 
algorithm. We will call this the ground plane baseline (even though it is not really a ground 
plane). The ground plane baseline is stronger than the baseline used in ( |Ashutosh Saxena 
and Ng|2007b"l ). 



We performed K-fold cross-validation using hard EM stereo pair training and view 
prediction testing. For each fold, 80% of the data was used for training and 20% was 
used for holdout testing. For testing, monocular inference was performed using loopy 



BP to minimize the energy function described in ( |5.2[ ). Some of the resulting monocular 
disparity maps for holdout test stereo pairs are shown in Figure |5.1[ demonstrating the 
ability of our monocular predictor to predict depth from different scenes. For the test 
pairs we also computed disparity maps using binocular only and combined monocular and 
binocular cues. The disparity prediction system combining both monocular and stereo cues 
appears to improve the accuracy of the disparity obtained from the initial stereo algorithm, 
although we did not have quantitative ground truth for the test pairs. 

The hard EM algorithm was run for four iterations of bootstrapping resulting in four 
weight vectors. Each monocular feature weight vector was tested by two measures on test 



data. First, we measured the view prediction error using formula ( |5.6 ) on the holdout stereo 



pairs. Second, we computed an error relative to ground truth on images from a second data 
set of 63 single images paired with laser range finder depth mapsjj We related disparity d 
to depth Z using Z = c/d where the constant c is determined by a fit to the entire data set 



3. We picked all the images with prefix 'img-stats' from the dataset described in ( Ashutosh Sax 



ena and Ng 2007b I. 
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Left image Stereo depthmaps Monocular depthmaps 

Figure 5.2: MDE results on some testing images 

(the same value of c is used for all test images). We reported RMS disparity error which is 
the square root of the average over pixels of (d — c/Z) . The results, averaged across all 
folds, are shown in Table |5.ll 



Table 5.1 shows that our algorithm converges quickly after only a few iterations. The 
prediction error improves for at least one iteration of training. This again suggests that 
the monocular cues improve the disparity estimation relative to binocular cues only. The 
table shows a dramatic improvement in RMS disparity error over the ground plane base- 
line. We obtained the average error in disparity estimate of 2 pixels, much better than the 
performance of the ground plane baseline (8.6 pixels). 

We also run our hard EM learning algorithm on the real-world road driving rectified 
stereo sequences acquired and provided by Daimler AG ( Liu and Klette|2007] ). The dataset 
is the composition of seven sequences, each of them presenting different traffic, road and 



94 
lighting conditions. Ground truth depth is not available. For this dataset, we trained our 
MDE model on six sequences and tested on the remaining sequence. As these are grayscale 
images, we used 15 image features at each scale. Again for testing, monocular inference 



was performed using loopy BP to minimize the energy function described in ( |5.2[ ). We 
demonstrated the resulting monocular disparity maps of a few frames in the testing se- 



quence in figure 5.2 The results show that the monocular model is able to provide good 
depth predictions from different scenes in general, although there were still image parts 
which were not generalized well. 

We have emphasized view prediction as a way of evaluating monocular depth esti- 
mation in the absence of any ground truth depth labels. We relate the learning problem 
using view prediction error to the problem of maximizing conditional likelihood in hard 
conditional EM learning. In Chapter [6| view prediction is employed again as an evalua- 
tion metrics for a structure and motion estimation framework. We have implemented our 
unsupervised CRF learning algorithm based on hard conditional EM for training MDE 
from stereo pairs only. The algorithm has shown to perform well by both view-prediction 
measures and by comparisons with ground truth labels, even though ground truth labels 
were not used in training. Furthermore, performance improved (modestly) over several 
iterations of learning which indicates that the monocular cues used in the EM process im- 
prove depth estimation when both monocular and binocular cues are used (as is reported in 
( |Ashutosh Saxena and Ng|2007a| )). In the second experiment, although the desired qualita- 
tive analysis could not be performed, we have demonstrated that MDE can also be used for 
Vision-based driving assistance. 



CHAPTER 6 
DEPTH AND MOTION ESTIMATION FROM STEREO 

SEQUENCES 

We introduce a unified framework for scene structure and motion estimation on road- 
driving stereo sequences (Trinh and McAllester 2010[ ). This framework is based on the 



slanted-plane scene model that has become widely popular in the stereo vision community. 
Our algorithm iteratively and alternately solves for scene structure and motion. Surface 
estimation is done using our own slanted-plane stereo algorithm. Motion estimation is 
achieved by solving a MRF labeling problem using Loopy Belief Propagation. We show 
that with some specific assumptions about the motion of the camera and the scene, the mo- 
tion estimation problem can be reduced to a ID search problem along the epipolar lines. 
We also propose a novel evaluation metrics, based on the notion of view prediction error 
previously introduced in Chapter [5} This metrics can be used to evaluate the performance 
of structure and motion estimation algorithms on stereo sequences without ground truth 
data. Experimental results on road-driving stereo sequences demonstrate that our algorithm 
successfully improves the view prediction error although it was not designed to directly op- 
timize this quantity. 

6.1 Introduction 

Multi-view video sequence is the richest form of image data for scene reconstruction. A 
subtype in this form that we consider very interesting are stereo video sequences. So far, 
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to our knowledge, this useful data type has hardly been investigated and exploited for the 
purpose of recovering scene structure and motion. Some of the very few research work that 



has discussed this issue was introduced by D. Min et al. (Min and Sohn 2006 1 and F. Huguet 
et al. ( |Huguet and Devernay|2007 ), in which the authors presented variational methods for 
scene flow estimation from calibrated stereo image sequences. This data type is interesting 
for two reasons. First, as opposed to a multi-view video sequence, a stereo sequence can 
easily be captured using only one moving calibrated binocular camera. Second, this data 
type provides us with enough constraints to conveniently compute location and motion of 
scene points simultaneously, since coupling dense stereo matching with motion estimation 
helps decrease the number of unknowns per image point. More specifically, for stereo 
sequences, we can formulate structure and motion estimation as an energy minimization 



problem, in which the model is either an extension of a stereo vision model ((Scharstein 



and SzeliskT||2002| )) to also handle scene point motion, or an extension of a Structure from 



Motion or optical flow model (( Simon Baker and Szeliski|2009 )) under a stereo setup. One 
example of using this latter formulation is the work presented in ( |Zhang and Kambhamettu 



2001). 



In this section, we follow the former approach. Based on our slanted-plane stereo 
model, we develop a new algorithm for structure and motion estimation on road-driving 
stereo sequences. This framework is based on the slanted-plane scene model that has be- 



come widely popular in the stereo vision community ( Scharstein et al. 2001 1. Our algorithm 
iteratively and alternately solves for scene structure and motion. Surface estimation is done 
using our own slanted-plane stereo algorithm. Motion estimation is achieved by solving a 
MRF labeling problem using Loopy Belief Propagation. We show that with some specific 
assumptions about the motion of the camera and the scene, the motion estimation problem 
can be reduced to a ID search problem along the epipolar lines. We also propose a novel 
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evaluation metrics, based on the notion of view prediction error. Again, it is also much eas- 
ier and more affordable to collect unlabeled stereo video data than to collect stereo video 
with associated ground truth (which usually requires using multiple sensors followed by a 
data fusion step). We argue that view prediction error can be used to evaluate the perfor- 
mance of structure and motion estimation algorithms on stereo sequences without ground 
truth data. Experimental results on road-driving stereo sequences support our argument by 
demonstrating that our algorithm successfully improve the view prediction error although 
it was not designed to directly optimize this quantity. We believe view prediction can be 
used as a quantitative performance measure on unlabeled multi-view datasets in a variety 
of applications. 

6.2 Proposed approach for Structure and Motion 

In this section, we present a detailed description of our unified formulation for scene struc- 
ture and motion. At each time step t in the video sequence we consider the stereo pair at 



time t and at time t + 1 (Figure 6.3 ). We only observe the three frames: II, In and it and 
the frame lj~ is unobserved. We want to estimate the 3-D location and motion for all pix- 
els in Ij^. To simplify the derivation of our framework, we use the following assumptions: 
the camera's viewing direction is the same direction of the Z axis, and all motion vectors 



in the scene are parallel to the Z axis. Later in section 6.3 we show that our approach still 
delivers good estimates even in sequences where these assumptions are violated, e.g: when 
there is small camera rotation, or rotation of moving objects. 

6.2.1 Connection between Disparity and Motion 

Our assumptions about the motion of the camera and the scene guarantee the following 
constraints: 
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The image location of the epipoles in all frames of the sequence is constant: As the 
epipole is simply the projection of the next camera center on the previous frame, this 
is obvious since the camera always moves along its viewing direction. 

From one frame to the next, a pixel translates along the epipolar line connecting that 
pixel and the epipole: since all 3D points in the scene move in the same direction 
with the Z axis, a 3D point always stay in the same 3D line which is parallel to the Z 
axis. Therefore the epipolar line, which is the projection of this 3D line on the image 



plane, stays constant, given that the epipole also stays constant. (Figure 6.1 ) 



Consider a 3D point P in the scene. Let R be the distance from the point to the Z axis. 
Let rt be the 2D distance from the projection of P to the epipole in the video frame at time 
t. We can derive the following mathematical relation between rt, rj+i, Vf and dt as follows 



(This relationship is also illustrated in Figure 6.2) 



n 

n+i 



R£ Rf Rd t 

z t fh/d t h 

Rf Rf 



z t+ i z t - Az t 

Rf Rd t 



fh/dt - Az t h(l - dtAzt/fh) 



n+i i i 



n 1 - dtAzt/fh 1 - dtvt 
dt+i n+i 1 



(6.1) 



(6.2) 



dt n 1 - d t v t 

where: / is the focal length, zt is the depth of P at time t (the actual distance of P w.r.t 

the camera), h is the stereo baseline, dt is the disparity of the projected pixel of P, vt is 
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the velocity value of P. Note that in the derivation above we make use of the following 
relationship between dt and zf. zt = fh/dp 

6.2.2 3-frame Model for Depth and Motion estimation 

At each time step t in the video sequence we consider the stereo pair at time t and the 



left frame at time t + 1. The frames are labeled as in Figure 6.3 By equation (6.1 ), we 
can see that if we know dt and vt, then we can compute rt+i, which means solving for 
the correspondence between /£ and I^~ . In other words, we converted the 2D optical 
flow field problem to the stereo disparity estimation problem and a ID velocity assignment 
problem. We define the following energy functions: 



E = E Stereo + E S + E M ( 6 - 3 ) 



Eg = ^min(r u ,A u \v P - vq\) (6.4) 

P,Q 



E M = ^2^2 X k 
P k 



4(p) 



fct+i 



(6.5) 



y - <f>k (p + Q(j>,dp,v p )) 



where Eg tereo is the function in (3.2), Q(p, d p , v p ) is the function mapping a pixel in 
l l L to a pixel in I^ , and 0*(p) is the k-dimensional feature vector corresponding to pixel p 
in frame at time t. The objective is to find the optimal co-assignment for depth and velocity 



that minimize the global energy function in (6.3 ). 



(d*,v*) = argmin(i^) 
(d,v) 
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The optimization of Eg tereo is introduced in section [4~T 

Minimizing the sum of the other two energy terms: E^j ot i on = E% + E V M can be 
considered another MRF labeling problem: we want to assign a velocity value to each su- 
perpixel such that Ey[ n on is minimized. The first term Eg penalizes the difference in 
velocity between two adjacent superpixels, the second term E V M is the data cost of assign- 
ing a velocity v p to pixel p, taking into account the image data agreement between frame 
l£ and frame J^ + . We solve this MRF labeling problem by Loopy Belief Propagation. 
We used the max-product BP algorithm with conceptually parallel updates. In our actual 
implementation however, they were performed sequentially. 

Finally, once we have obtained an initial estimate of vt, we can fix Vf and optimize 
the function Ej, rame = Egtereo + E V M . Note that this function has 1 more data term 
compared to the previous stereo function as we can now incorporate data from I^~ . We 
construct an iterative algorithm alternately optimizing the disparity dt and the velocity vt 
for frame 1^. Here is the outline of our algorithm: 

1 . Compute the epipole at it using a version of the 8-point algorithm. 

2. Estimate the disparity map dt of I l L using our stereo algorithm. 

3. Use SIFT feature matching to obtain a set of initial velocity values: Each pair of 
matched SIFT features give us 1 initial velocity. 

4. Estimate vt given df. Run Loopy BP to optimize Ej\^ ot ^ on and assign a velocity value 
to each superpixel in /£. 

5. Fixing vt, reestimate dt by optimizing EJ ereo . 

6. Repeat from Step 4. 
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6.3 Experimental Results 

6.3.1 Evaluation Metrics 

We use view prediction as evaluation metrics for our experiments here. The idea of using 



view prediction error has been investigated in Section 5.2 in which we described stereo 
pair view prediction. View prediction is a general notion that has been motivated by the 
two-view learning approach in machine learning (Blum and Mitcheilj |1998 ). In two-view 



learning, the goal is to predict the second view given the first. In our context the goal of 
view prediction is to predict the unobserved data y given the observed data x and the latent 
variables w.We can express a probabilistic interpretation of view prediction as follows: 



* 
w 



argmin \^ In 1/P w (y\x) 



where P w (y\x) is the conditional probability of the unobserved data given the observed 
data parametrized by w. As our algorithm never observed /p" , in our problem setting 
we might define P w {y\x) as the P{I^ 11^,1^,1^ ,dt,vt). The idea is that the more 
accurate we estimate (dt, Vf), the more accurate we can predict I p . 

Given the estimated depth and velocity dt, Vf, we can compute the prediction In as 



follows: 



i\ dt 4 l 4+ 1 ^ IJ+ 1 



where c^+i is computed using (6.2). An example of a predicted fourth view is shown 



in Figure 6.4 



The view prediction error is defined as the pixel RMS between I R + and I R + : 
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Err(I t + L ,I t + i ) = \I ^P= l[R ^ R±1L. (6 . 6) 



r m f t + u J SLi(T(p)-T(p)) ! 

l R i'r ) ~ V N 



We can see a clear analogy between equation (6.6) and equation (5.6) from Section 



5.2 



The difference is that the computation of the prediction In here is more complicated, 



involving both the estimated disparity and velocity. In the absence of ground truth data, it is 
natural to use view prediction error as a useful tool for quantitative analysis, i.e. to measure 
how good the estimation of latent variables is. In the specific setting of our problem, the 
latent variables w* = (d*,v*). 

6.3.2 Results on road-driving stereo sequences 

We tested our algorithm on 7 gray scale road-driving stereo sequences provided by Daimler 
AG. Each sequence contains from 250 to 300 rectified, bias gain corrected stereo pairs, 
taken under different light condition and road setting. Ground truth depth and motion is 
not available for these datasets. The algorithm estimates the dense map of disparity and 
velocity value for each left frame in the sequence. Figure [63 demonstrates our results on 
several frames in a sequence. 



Table 6.1 shows the average view prediction error ( 6.6 ) over seven video sequences after 
3 iterations of the algorithm. The results shows that the algorithm successfully improved 
the estimation over time. The improvement stops after three iterations. Note that the pixel 
intensity values are scaled to be in the range of [0 — 1]. 
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Average View Prediction Error 


Iter 1 


0.0692 


Iter 2 


0.0622 


Iter 3 


0.0621 



Table 6.1: Average view prediction error (in pixel intensity value) on 7 Daimler road- 
driving stereo sequences after each iteration of the algorithm. 

6.4 Unsupervised Learning of the Three-frame Model 

In Chapter |4j we demonstrated our algorithm for unsupervised learning of a highly pa- 
rameterized stereo vision model involving the shape from texture cues - training the model 
parameters from unlabeled stereo pair training data. Our approach to unsupervised learning 
is based on maximizing conditional likelihood. 

Throughout Chapter [4j we used stereo vision as a simple setting to investigate un- 
supervised learning. However we also argued that our approach to unsupervised learning 
based on maximizing conditional likelihood can be generalized to other, maybe much more 
sophisticated models. 

In this chapter, we formulate our unsupervised learning algorithm specifically for the 



model described in Section 6.2 We rewrite the 3-frame depth and motion estimation model 
here: 
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E = E Stereo + E S + E M ( 6 - 7 ) 



E S = ^min(r v ,A v \v P - vq\) (6.8) 

6.4.1 Unsupervised training using Hard Conditional EM/Bootstrap 

training 

In our three-frame model, we denote x to be the frame it, y = {y^\y^ '), where y^ > 
is the frame I f R , y^> is the frame Iw~ , and z = (d, v) to be the latent variable, i.e the 
assignment of disparity and velocity to all pixels in x. 



We are interested in training the parameter vector j3 = (/3^, (3 V ) of the model in (6.7 ) so 
as to optimize the conditional probability of y given x. 



N 
(3* = argmax ^ ^PpiVi^i) ( 6 - 10 ) 

N 

= argmax ^ mP /?(l/l >y\ \ x i) ( 6 - n ) 

I 3 i=l 



Hard conditional EM locally optimizes the following version of (6.10) 



N 



ff = argmax \. max m -P/3(?/i> z \ x i) 



i=l 



N 



argmax N. max m -P/3(^ iVi id-> v \ x i) 



Specifically, hard conditional EM iterates the following two updates. 
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(6.12) 



(6.13) 



zi := axgmaxP^(yi,z\xi 



(di,Vi 



P 



D I (!) ( 2 ) A I 

argmax Pq\v\ ' ,y\ \d,v\xi 

(d,v) 

N 
argmax V"lnP«(^,^|^) 



N 



Wd, Pi 



.(1) J 2 ) 



argmax $^ lnP (/3 d A>)^ '^ , d h v iWj 



(6.14) 
(6.15) 

(6.16) 

(6.17) 



Equation (6.14) is the hard E step and equation ( 6.16| ) is the hard M step. In the case 
of our three-frame model, the hard E step is implemented using the inference algorithm 



described in Section 6.2 This iterative algorithm finds the optimal co-assignment for depth 
d and velocity v that minimize the global energy function in ( |6.3[ ). 

To implement the hard M step, we first rewrite the probability in the right hand side of 



(6.16) as follows. 



P(/3 d A)(y {1 \y {2 \d,v\x) = P^{y^ : d\x)P^ 2 \v\d : x) 



(6.18) 
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The first term, Pa (yd),d\x) can be further factorized as demonstrated in (4.30), (4.31 ) 



and (4.32). 



The second term, Pa (y^>,v\d, x) can also be factorized in a very similar way. 



P Mv {y {2 \v\d,x) = Pa(yW\v,d,x)Pp v (v\d,x) 



(6.19) 



PB v (y (2) \v : d, 



X) 



e (- E M&y,v,P y ) ) 

Zy(x,V,(3y) 



(6.20) 



Zy( X ,V,Py) = ^(-^(WA)) 



Pp v (v\d,x) 



e (-E v s (x,v,P v ) ) 
Z v (x, p v ) 



(6.21) 



Z v (x,/3 V ) = ^ 



e (-E%(x,v,/3v) ) 



Given this factorization, the hard M step in equation (6.16) then can be implemented 
the same way as described in equations ( |4.33[ ), ( |4.33[ ) and ( |6.22[ ). 

Specifically, let E{(v) abbreviate Eg(x, v,(3 v ),the contrastive divergence update equa- 
tion for j3 v can be expressed as follows. 



N 
V /3, L = E (Ev-MvM i V Pv E i( v )} - V^Eiivi)) (6.22) 



i=\ 



A similar equation holds for /3 y and E v M (x, y, v, /3 y ). 



107 
6.4.2 Unsupervised training using fourth view prediction 



Here we are interested in training the parameter vector of the model in (6.7 ) so as to min 



imize the fourth view prediction error (6.6). First, assume that we have a way to compute 



the gradient of the view prediction error function (6.6) with respect to the latent variables 
d, v. Using contrastive divergence can be one of the possibilities. We call this gradient G. 
Second, we can also assume that the gradient of the latent variables d and v with respect to 
the model parameters (3 can also be computed. Let g be this gradient. 

By applying the chain rule, the gradient of the view prediction error function with 
respect to /3 can be computed as follows. 

Q0 R = G 9 (6-23) 

We can then use gradient descent to find the optimal j3* optimizing the fourth view 



prediction error (6.6 ). 



6.5 Conclusion 

In this chapter we emphasized view prediction as a way of evaluating scene structure and 
motion estimation from stereo video data in the absence of any ground truth labels. We 
introduced a new algorithm for structure and motion estimation on road-driving stereo se- 
quences. Based on specific assumptions about the motion of the camera and the scene, we 
can reduced the 2D optical flow problem to a ID velocity value problem. Our algorithm it- 
eratively and alternately solve for structure and motion. Scene structure estimation is done 
using our own plane-based stereo algorithm. Velocity estimation is completed by solving 
a MRF labeling problem using Loopy BR Experiments on road-driving stereo sequences 
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showed encouraging results, even with video sequences where the scene and camera mo- 
tion do not fully comply with our assumptions. 

Performance analysis was done using our novel evaluation metrics based on the no- 
tion of view prediction error. We argue that this evaluation metrics is quite appropriate 
for algorithms working with stereo sequences as well as multiple view image data, when 
ground truth data is not available. Our experimental results support this argument. We used 
hand-tuned parameters for our model in this chapter. Ideally, these parameters should be 
estimated by an automated method, usually through learning using a labeled training data 
set. With view prediction error, we believe the problem we solve in this chapter can be 
another setting where we can apply the unsupervised learning approach based on maximiz- 
ing conditional likelihood. In other words, the model parameters can be learned using only 
unlabeled stereo video data. 

Since we could relate the learning problem using view prediction error to the problem 
of maximizing conditional likelihood in hard conditional EM learning, as demonstrated 
in Chapter [5} in addition to the fact that the problem we address here is also modeled as 
CPvFs, there is obviously the possibility of using this problem as another setting where we 
can apply our unsupervised CRF learning approach. It is hoped that a general notion of 
view prediction will eventually facilitate unsupervised learning in a variety of cases with 
mutual information between views. 
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Figure 6.1: As the camera moves forward, each pixel moves along its corresponding epipo- 
lar line. The distance it moves depends on both their depth and velocity. Since the truck on 
the left has highest velocity w.r.t the camera, P\ moves the longest distance. On the other 
hand, the truck on the right has almost zero velocity w.r.t the camera, hence P<i almost 
stands still. P3 reflects the motion of the camera w.r.t the road (the scene background). 
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Figure 6.2: The geometric interpretation of equation (6.1 ). 
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Figure 6.3: 
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Figure 6.4: Top: original fourth frame / jt . Bottom: predicted fourth frame In . The 
black pixels are missing values caused by bilinear interpolation. These pixels are excluded 
when computing the view prediction error. 
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Figure 6.5: Results on a Daimler road driving sequence. 



CHAPTER 7 
CONCLUSION 

In this thesis, we demonstrated machine learning techniques to build computer vision sys- 
tems 3-D scene geometry recovery from images. Scene structure can be recovered from 
different types of input image data. In this thesis we focused on three data types: two-view 
stereo pairs, single image, and stereo video sequences. Building these systems requires 
algorithms for doing inference as well as learning the parameters of MRFs or CRFs. 

For inference in systems with continuous valued variables, or discrete-valued variables 
with very large domains, it is impossible to directly use a standard discrete MRF inference 
techniques such as Loopy BP, graph cuts or tree-reweighted message passing. For such 
systems, we proposed to use a generic Particle-based Belief Propagation (PBP) algorithm 
closely related to previous work, which we then formulated specifically for our MRF label- 
ing problems. Although we only described a specific use of this generic PBP algorithm, we 
believe it can be used as an approximate inference scheme for a wide variety of problems 
that can be formulated by a probabilistic graphical model, especially those containing many 
random variables with very large or continuous domains. Our approach creates a set of par- 
ticles for each variable, representing samples from the posterior marginal. The algorithm 
then continues to improve the current marginal estimate by constructing a new sampling 
distribution and draw new sets of particles. Its shown by experiments that the algorithm 
is consistent, i.e. approaches the true values of the message and belief functions with fi- 
nite samples. Besides accuracy, PBP algorithm also provides good estimate for properties 
of the distribution, and a representation of state uncertainty. Although an adaptive choice 
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for the sampling distribution, and such iterative resampling processes require further work 
to analyze, our results seem to support the notion of sampling from the current marginal 
estimates themselves, whether from fitted distributions or via a series of MCMC steps. 

Learning in MRF involves estimating the optimal parameters of the energy model. For 
learning, unlike previous work, we trained our system without using ground-truth labeled 
data. We introduce the unsupervised CRF learning algorithm, based on maximizing con- 
ditional likelihood using hard conditional EM. We demonstrated the application of our 
unsupervised CRF learning algorithm for different scene geometry problems 

In unsupervised learning one usually formulates a parameterized probability model and 
seeks parameter values maximizing the likelihood of the unlabeled training data. Specif- 
ically for stereo vision, the model defines the conditional probability of the right image 
given the left. We introduced a slanted-plane stereo vision model in which we used a fixed 
over- segmentation to segment the left image into coherent regions called superpixels. We 
then assign a disparity plane for each superpixel. We formulated the problem of inferring 
plane parameters as a MRF labeling problem, which can be solved by an energy mini- 
mization method. The MRF is a graphical model in which superpixels define nodes and 
the adjacency between superpixels define edges. Our stereo energy function balances be- 
tween a data matching term and a smoothness term. We then used our learning algorithm 
to train this highly parameterized stereo vision model involving the shape from texture 
cues. We showed that this unsupervised learning algorithm implicitly trains shape from 
texture monocular surface orientation cues. We demonstrated that training monocular cues 
from stereo pair data improved stereo depth estimation. Our stereo model with texture 
cues, only by unsupervised training, outperformed the results in related work on the same 
stereo dataset. Stereo vision provides perhaps the simplest setting in which to study un- 
supervised learning. We have formulated an approach to unsupervised learning based on 
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maximizing conditional likelihood and demonstrated its use for unsupervised learning of 
stereo depth with monocular depth cues. Ultimately we are interested in learning highly pa- 
rameterized sophisticated models including, perhaps, models of surface types, shape from 
shading, albedo smoothness priors, lighting smoothness priors, and even object pose mod- 
els. We believe that unsupervised learning based on maximizing conditional likelihood can 
be scaled to much more sophisticated models than those demonstrated in this paper. 

We also applied our unsupervised CRF learning approach to the problem of monocular 
depth estimation (MDE). We learned the MDE model using stereo pairs without ground 
truth depth maps. While most of the performance of the monocular predictor is achieved 
by training on the initial stereo depth estimates, a modest improvement is seen for both view 
prediction and for ground truth prediction at later EM iterations. The MDE performance is 
similar to that reported in previous work on the same dataset despite the absence of ground 
truth in learning. 

In this thesis, we also addressed the use of an interesting image data type - stereo 
video sequences. Specifically, for stereo sequences, we can formulate structure and motion 
estimation as an energy minimization problem, in which the model is either an extension of 
a dense stereo vision model to also handle scene point motion. In this thesis, based on the 
constructed slanted-plane stereo model, we introduced a new algorithm for structure and 
motion estimation on road-driving stereo sequences. Based on specific assumptions about 
the motion of the camera and the scene, we can reduce the 2D optical flow problem to a ID 
velocity value problem. Our algorithm iteratively and alternately solve for structure and 
motion. Surface estimation is done using our own slanted-plane stereo algorithm. Velocity 
estimation is achieved by solving a MRF labeling problem using Loopy BR Performance 
analysis was done using our novel evaluation metrics based on the notion of view prediction 
error. Experiments on road-driving stereo sequences showed encouraging results, even 
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with video sequences where the scene and camera motion do not fully comply with our 
assumptions. We emphasized view prediction as a way of evaluating scene structure and 
motion estimation from stereo video data in the absence of any ground truth labels. We 
argue that this evaluation metrics is quite appropriate for algorithms working with stereo 
sequences as well as multiple view image data, when ground truth data is not available. 
Our experimental results support this argument. We used hand-tuned parameters for our 
model in this paper. Ideally, these parameters should be estimated by an automated method, 
usually through learning using a labeled training data set. Our thesis addressed a theoretical 
discussion regarding this issue, although we did not actually implement the idea. With view 
prediction error, we believe the scene structure and motion estimation from stereo video 
problem can be another setting where we can apply our unsupervised learning approach 
based on maximizing conditional likelihood. In other words, the model parameters can be 
learned using only unlabeled stereo video data. 
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